RNA-Seq Data Visualization & Analysis with R
Matt Maxwell
April 15, 2024
RNA-seq analysis pipeline
The RNA-seq analysis pipeline sections colored in gray is not covered here but were performed using HOMER on the Salk Institute IGC high performance computing cluster. The pipeline sections highlighted in pink are performed in Rstudio and are covered in depth here including the R code and resulting outputs. RNA-seq analysis and data visualization have been critical discovery tools in my doctoral work seeking to define mechanisms by which mutations in the ARID1A tumor suppressor can influence anti-tumor immunity. For example, via GSEA analyses I found that Type I Interferon gene signature activation is commonly observed in cancer cell lines following ARID1A deletion.
## Warning: package 'htmlwidgets' was built under R version 4.3.2Import data & load packages
Data analyzed here is from an RNA-seq experiment comparing the
transcriptomes of MC38 mouse colon cancer cells deficient in the ARID1A
tumor suppressor gene (sgArid1a) relative to CRISPR control
MC38s (sgScramble). The MC38_edgeR_df dataframe from this
experiment includes important data such as normalized mRNA transcript
counts, Log2fc, and adjusted p values across the transcriptome outputted
from edgeR.
#To set working directory on windows machine, replace backslashes with double backslashes for R to accept our file path
wd <- r"(C:\Users\mattm\OneDrive\Desktop\GitHub_projects\RNA_Seq_Data_Visualization_with_R)"
#set working directory
setwd(wd)
#set this as our working directory for all R code chunks in this tutorial.
#IGNORE if you're not using R markdown file
knitr::opts_chunk$set(root.dir = wd)
#Set preference for no warnings to be given when running code chunks in R markdown
#IGNORE if you're not using R markdown file
knitr::opts_chunk$set(warning = FALSE, message = FALSE)
#load packages with package manager package called pacman
if (!require("pacman")) install.packages("pacman")
pacman::p_load(here,
tidyverse, #dyplyr and ggplot utilities
scales, # Transform axis scales
ggrepel, # Optimize plot label separation
ComplexHeatmap, #Awesome heatmap package
grid, #plot heatmaps side by side
viridis, #Package for some of my favorite color palettes like Magma
RColorBrewer,#Package for some of my favorite color palettes like RdBu
forcats, #Package for working with factors
msigdbr, #Molecular signatures database package for GSEA
clusterProfiler, #Package for GSEA dotplot
euler) #Package for making venn diagramsData wrangling for volcano plot
To make a volcano plot, we need three pieces of data: 1) gene name, 2) Log2 fold change, and 3) adjusted p values. Here, I demonstrate how to grab these three pieces of data and use them make a volcano plot that highlights genes whose expression significantly changes as well as to specifically highlight some individual genes I think are interesting.
#Import data
MC38_edgeR_df <- read.delim("sgArid1a_MC38_edgeR.txt", header = TRUE, sep = "\t")
#Rename our columns of interest using the rename function in dyplyr from tidyverse package
MC38_edgeR_df <- rename(MC38_edgeR_df, Gene = Annotation.Divergence, Log2fc = sgArid1a.vs..sgScramble.Log2.Fold.Change, adj_p_value =sgArid1a.vs..sgScramble.adj..p.value)
#Grab the gene name values preceding "|" edgeR from Gene column
MC38_edgeR_df$Gene <- sub("\\|.*", "", MC38_edgeR_df$Gene)
#grab the columns we need for volcano plot
volcano_df <- MC38_edgeR_df[, c("Gene", "Log2fc", "adj_p_value")]
# Replace NA values in Log2FC column with Zero.
volcano_df$Log2fc[is.na(volcano_df$Log2fc)] <- 0
#Define the parameters you want to classify a gene as upregulated, donwregulated, or not significantly changed (NS)
volcano_df <- volcano_df %>%
mutate(gene_type = case_when(Log2fc >= .585 & adj_p_value <= 0.05 ~ "Upregulated",
Log2fc <= -.585 & adj_p_value <= 0.05 ~ "Downregulated",
TRUE ~ "NS"))
#Count the number of genes in the three classes we defined
volcano_df %>%
count(gene_type)## gene_type n
## 1 Downregulated 218
## 2 NS 24135
## 3 Upregulated 184#Specify colors, sizes, and transparancy values associated with the three classes
cols <- c("Upregulated" = "#ffad73", "Downregulated" = "#26b3ff", "NS" = "grey")
sizes <- c("Upregulated" = 3, "Downregulated" = 3, "NS" = 1)
alphas <- c("Upregulated" = 1, "Downregulated" = 1, "NS" = 0.5)
# Define genes to highlight on the volcano plot
signif_genes <- volcano_df %>%
filter(Gene %in% c("Tap1", "Ifi27l2a", "Stat1", "Cxcl10", "Ptgs1", "Cd47"))
Up_genes <- volcano_df %>%
filter(Gene %in% c("Tap1", "Ifi27l2a", "Cxcl10", "Stat1"))
Down_genes <- volcano_df %>%
filter(Gene %in% c("Ptgs1", "Cd47"))Volcano plot
Here, we use ggplot2 from the Tidyverse R package to make the volcano plot because ggplot allows more room for customization when making a volcano plot compared to some other R packages. Volcano plots such as this are a great way to visualize pair wise comparisons of the number of genes changing, the magnitude of change, and statistical confidence associated with those changes. We’ll also be sure to color code genes which we define as differentially expressed genes (DEGs) based on their log2fc and log10p values.
#Customize what we want to show in our the volcano plot
volcano_df %>%
arrange(factor(gene_type, levels = c("NS", "Upregulated", "Downregulated"))) %>% #arrange gene_type in desired order for plotting!
ggplot(aes(x = Log2fc, y = -log10(adj_p_value))) +
geom_point(aes(colour = gene_type),
alpha = .75,
shape = 16,
size = 4) +
geom_point(data = Up_genes,
shape = 21,
size = 4,
fill = "firebrick",
colour = "black") +
geom_point(data = Down_genes,
shape = 21,
size = 4,
fill = "steelblue",
colour = "black") +
geom_hline(yintercept = -log10(0.05),
linetype = "dashed") +
geom_vline(xintercept = c(log2(.75), log2(1.5)),
linetype = "dashed") +
geom_label_repel(data = signif_genes,
aes(label = Gene),
force = 1,
nudge_y = 2) +
scale_colour_manual(values = cols) +
scale_x_continuous(breaks = c(seq(-6, 8, 2)),
limits = c(-6, 8)) +
labs(title = "MC38 Volcano Plot",
x = "Log2FC \n(sgArid1a/sgScramble)",
y = "-Log10(adj P-value)",
colour = "Gene type") +
theme_bw() + # Select theme with a white background
theme(panel.border = element_rect(colour = "black", fill = NA, linewidth = 0.5),
panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
axis.title.x = element_text(size = 11),
axis.title.y = element_text(size = 11),
plot.title = element_text(size = 14, face = "bold", hjust =.425))
For presenting data, it is helpful to highlight genes of interest on our
volcano plot such as the xix genes I’ve highlighted on this volcano
plot. For context, the two downregulated DEGs I highlighted are
immunosuppressive and four upregulated DEGs are associated with an
inflammatory interferon response which could contribute to our
concurrent observation of anti-tumor immunity in sgArid1a MC38
tumors compared to sgScramble MC38 tumors.
Data wrangling for MA plot
While a volcano plot visualizes the Log2FC and significance of DEGs, it tells you nothing about the expression level of these DEGs. An MA plot can fill this gap by plotting RNA-seq data normalized counts (TPM) vs Log2FC. In this code chunk, we join the mean TPM values of genes together with their Log2FC for plotting.
#read in TPM
tpm <- read.delim("MC38_sgArid1a_TPM_.txt", header = T)
#Clean gene names
tpm$Annotation.Divergence <- sub("\\|.*", "", tpm$Annotation.Divergence)
#Rename gene column
tpm <- tpm %>%
rename(Gene = Annotation.Divergence)
#Calculate mean TPM
tpm <- tpm %>%
rowwise() %>%
mutate(
sgNT_mean = mean(c_across(starts_with("sgNT")), na.rm = TRUE),
sgArid1a_mean = mean(c_across(starts_with("sgArid1a")), na.rm = TRUE),
mean_TPM = mean(c(sgNT_mean, sgArid1a_mean), na.rm = TRUE) # Calculate the mean of sgNT_mean and sgArid1a_mean
) %>%
ungroup()
#select gene and tpm means for sgNT and sgArid1a
tpm_means <- tpm %>%
select(Gene, sgNT_mean, sgArid1a_mean, mean_TPM)
#Merge dataframes to include mean_TPM column
MA_df <- left_join(tpm_means, volcano_df, by = "Gene")
Up_genes <- left_join(Up_genes, tpm_means, by = "Gene")
Down_genes <- left_join(Down_genes, tpm_means, by = "Gene")
signif_genes <- left_join(signif_genes, tpm_means, by = "Gene")MA plot
An immediate insight from the MA plot below is that the gene with the largest Log2FC of sgArid1a upregulated genes, Ifi27l2a is expressed at a much lower level than the other highlighted upregulated DEGs. Further inspection of this data reveals this gene is turned off in sgScramble cells but highly indued in sgArid1a cells, an interesting finding that’s facilitated via MA plot.
#MA Plot
MA_df %>%
arrange(factor(gene_type, levels = c("NS", "Upregulated", "Downregulated"))) %>% #arrange gene_type in desired order for plotting!
ggplot(aes( x = log2(mean_TPM),
y = Log2fc)) +
geom_point(aes(colour = gene_type),
alpha = 1,
shape = 16,
size = 3) +
geom_point(data = Up_genes,
shape = 21,
size = 4,
fill = "firebrick",
colour = "black") +
geom_point(data = Down_genes,
shape = 21,
size = 4,
fill = "steelblue",
colour = "black") +
geom_hline(yintercept = 0,
linetype = "dashed") +
geom_label_repel(data = signif_genes,
aes(label = Gene),
force = 1,
nudge_y = 2) +
scale_colour_manual(values = cols) +
labs(title = "MC38 MA Plot",
x = "Log2(Mean TPM)",
y = "Log2FC \n(sgArid1a/sgScramble)",
colour = "Gene type") +
theme_bw() + # Select theme with a white background
theme(panel.border = element_rect(colour = "black", fill = NA, linewidth = 0.5),
panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
axis.title.x = element_text(size = 11),
axis.title.y = element_text(size = 11),
plot.title = element_text(size = 14, face = "bold", hjust =.45))
Data wrangling for Heatmap
To make a heatmap, we need to 1) grab the normalized counts values for each sample from our edgeR dataframe (MC38_edgeR_df), 2) filter for genes we want to visualize such as genes whose expression significantly changes, and 3) convert this subset of genes from a dataframe to a matrix to use as input for the ComplexHeatmap. In addition, we will also utilize some code to highlight genes of interest (same genes of interest on volcano plot) on the heatmaps we make.
#Rename the normalized counts columns we'll need for our heatmap by referencing their column number
MC38_edgeR_df <- rename(MC38_edgeR_df, sgScr_1 = 9, sgScr_2 = 10, sgArid1a_1 = 11, sgArid1a_2 = 12)
#Filter for significantly up or downregulated genes
#We only want to plot significantly changed genes on this heatmap
MC38_edgeR_df2 <- MC38_edgeR_df %>%
filter(Log2fc >= .585 |Log2fc <= -.585, adj_p_value <= .05)
#Grab gene name and edgeR normalized counts for each group
heatmap_df <- MC38_edgeR_df2[, c(8:12)]
#Make dataframe without gene column, the first column of the df
heat_mat <- heatmap_df[,-1]
#Assign gene names from 'heatdata' as row names to new dataframe object 'mat'
rownames(heat_mat) <- heatmap_df[,1]
#Convert dataframe to a matrix, can use data.matrix() or as.matrix()
heat_mat <- data.matrix(heat_mat)
#Generate gene Z scores and transpose matrix
heat_mat <- t(scale(t(heat_mat)))
#List of interesting genes we'd like to label on heatmap
heat_anotation <- signif_genes
#add a column for row # in heatmap_df
heatmap_df$rownumber = 1:nrow(heatmap_df)
#Add row number values from heatdata to heat_anotation
heat_anotation <- merge(heat_anotation,heatmap_df, all = F)
#Print row numbers for heat_anotation from heatdata
cat(heat_anotation$rownumber,sep=",")## 229,234,303,352,295,16x <- heat_anotation$Gene
x_list <- as.list(x, ",")
print(x_list)## [[1]]
## [1] "Cd47"
##
## [[2]]
## [1] "Cxcl10"
##
## [[3]]
## [1] "Ifi27l2a"
##
## [[4]]
## [1] "Ptgs1"
##
## [[5]]
## [1] "Stat1"
##
## [[6]]
## [1] "Tap1"#Create object with row name locations of genes from heatdata that you'd like to annotate on heatmap
genelabels = rowAnnotation(foo = anno_mark(at = c(229,234,303,352,295,16),
labels = x_list))Heatmaps
To make the heatmap from our data matrix, we simply need to use the Heatmap function from the ComplexHeatmap package. However, I also want to make some additional asethetic modifications to my heatmap such as giving each genotype their own color coded top label and specifying that the heatmap utilize the “RdBu” color palette from the RColorBrewer R package in the first heatmap or the “Magma” palette from the Viridis package in the second heatmap.
#Define a list of our sample genotypes as a factor of either sgScramble or sgArid1a
fa = factor(c("sgScramble", "sgScramble", "sgArid1a", "sgArid1a"),
levels = c("sgScramble", "sgArid1a"))
#Define Heatmap annotation color bars for genotype designation
ha = HeatmapAnnotation(Genotype = fa, height = unit(0.5, "cm"),
col = list("Genotype"=c("sgScramble"="black","sgArid1a"="hotpink")))
#Passing my fav color palette RdBu to mypalette object
mypalette = rev((brewer.pal(n=9, ("RdBu"))))
ht1 <- Heatmap(heat_mat,
col = mypalette,
top_annotation = ha, heatmap_width = unit(5 , "cm"), heatmap_height = unit(15, "cm"),
border = TRUE, cluster_rows = T, right_annotation = genelabels,
show_column_dend = F, show_row_dend = F, show_row_names = F, column_names_side = "top",
column_names_rot = 2, show_column_names = F, cluster_columns = F,
use_raster=T, row_km = 1, raster_quality=5,
column_names_gp = gpar(fontsize=8),
heatmap_legend_param = list(height = unit(8, "cm"), direction = "horizontal",
title = "Expression Z Score", border = "black",
title_position = "topleft"))
ht2 <- Heatmap(heat_mat,
col = viridisLite::magma(n = 100),
#col = mypalette,
top_annotation = ha, heatmap_width = unit(5 , "cm"), heatmap_height = unit(15, "cm"),
border = TRUE, cluster_rows = T, right_annotation = genelabels,
show_column_dend = F, show_row_dend = F, show_row_names = F, column_names_side = "top",
column_names_rot = 2, show_column_names = F, cluster_columns = F,
use_raster=T, row_km = 1, raster_quality=5,
column_names_gp = gpar(fontsize=8),
heatmap_legend_param = list(height = unit(8, "cm"), direction = "horizontal",
title = "Expression Z Score", border = "black",
title_position = "topleft"))
ht_list <- ht1 + ht2
draw(ht_list, heatmap_legend_side = "bottom")
Heatmaps are useful for visualizing the relative relationships between gene expression values between samples where rows are individual genes and columns are samples. You can also communicate how many genes are upregulated versus downregulated in relative terms with heatmaps. Here, I’ve highlighted the same six genes I highlighted in the volcano plot that are implicated inthe biological process I’m studying, anti-tumor immunity.
Data wrangling for GSEA
To provide an informative visualization of which biological processes are changing in an experimental condition, we’ll use gene set enrichment analysis (GSEA) to identify which gene sets are enriched in in sgArid1a or sgScramble MC38 cells. Note: an important distinction between GSEA and pathway overepresentation analyses (ORA) is that GSEA is performed on all genes while ORA is performed on a filtered list of genes such as DEGs. This is powerful because GSEA can detect subtle but significant enrichments in gene sets across the transcriptome that may have been missed by filtering DEGs before pathway ORA analysis. To begin, we grab gene names, Log2FC, and adjusted p values for all our genes. We then assign genes a GSEA score and sort them based on their GSEA score to make a ranked list in descending order. We then perform GSEA analysis on the ranked list using the ClusterProfiler R package and plot the results using ggplot.
#Grab gene name, Log2fc, and adj_p_value
GSEA_df <- volcano_df[,1:3]
#Add a column for GSEA values
GSEA_df$GSEA_score <- -log10(GSEA_df$adj_p_value) * sign(GSEA_df$Log2fc)
#Create ranklist object with gene name and GSEA value
rankList <- GSEA_df[,4]
names(rankList) <- as.character(GSEA_df[,1])
#sort the genes in the list according to highest to lowest GSEA value
rankList <- sort(rankList,decreasing=TRUE)
#Lets see the top genes in the list
head(rankList)## Slc5a3 Ifi27 Cotl1 Helz2 Lgals3bp Oasl2
## 65.46435 61.78099 53.29668 51.30068 49.41965 41.24129#Use msigdbr package to read in the gene sets for GSEA analysis
#msigdbr gives access to many different gene sets
#For example, instead of Hallmarks, could import GO gene sets by specificying category ="C5"
Hallmark <- msigdbr(species = "Mus musculus", category = "H")
Hallmark_v2 = Hallmark %>% dplyr::select(gs_name, gene_symbol) %>% as.data.frame()
#Run GSEA analysis of your gene list against Hallmarks
set.seed(8888)
Gsea_result <- GSEA(rankList, TERM2GENE=Hallmark_v2, verbose=FALSE,
pvalueCutoff = .8, pAdjustMethod = "BH", nPerm = 1000)
#Manipulate strings (eg., str_to_title function)
library(stringr)
## count the gene number for each enrichment
Gene_count<- Gsea_result@result %>% group_by(ID) %>%
summarise(count = sum(str_count(core_enrichment,"/")) + 1)
## merge with the original dataframe
GSEA_data_frame <- left_join(Gsea_result@result, Gene_count, by = "ID") %>% mutate(GeneRatio = count/setSize)
## for reordering the factor
library(forcats)
Data_activated1 <- GSEA_data_frame %>% filter(NES>0 | NES < 0)
#Take the word Hallmark out of gene sets for plot visualization purposes
Data_activated1$Description <- gsub("_", " ", Data_activated1$Description)
Data_activated1$Description <- str_to_title(Data_activated1$Description)
Data_activated1$Description <- gsub("Hallmark ", "", Data_activated1$Description)
#Arrange gene set df by NES score
Data_activated1 <- arrange(Data_activated1, desc(NES))
#Take the hightest and lowest NES value gene sets
upregulated_gs <- slice_head(Data_activated1, n=5)
downregulated_gs <- slice_tail(Data_activated1, n=5)
#Combine the upregulated_gs and downregulalted_gs objects to plot on dotplot
Dotplot_gs <- rbind(upregulated_gs, downregulated_gs)GSEA Dotplot
Plot GSEA dotplot with ggplot.
custom_dotplot <- ggplot(Dotplot_gs, aes(NES, fct_reorder(Description, NES))) +
geom_point(aes(size = GeneRatio, color = p.adjust)) +
theme_minimal(base_size = 15) +
#scale_colour_viridis_c(limits = c(5.0e-06, 2.5e-05). direction = 1, option = "viridis")+
ylab(NULL) +
#ggtitle("Hallmark enrichment", ) +
scale_size_continuous(range = c(3, 14)) + scale_colour_gradient(low="red", high = "blue") + coord_cartesian(clip = "off") +
labs(color="Adjusted p-value", size="Gene Set Ratio") +
theme(axis.text=element_text(size=10,color="black")) +
xlab("NES (sgArid1a/sgScramble)") +
ggtitle(substitute(paste(bold("GSEA Hallmarks Dotplot")))) +
theme(axis.title.x = element_text(size = 15, angle = 0, vjust = -1))
custom_dotplot
Here we can see the top enrichments in sgArid1a MC38 are
pro-inflammatory Interferon Responses, a central element of my thesis
studies which are known to spark anti-tumor immunity.
Proportional Venn Diagrams
We can also make venn diagrams to assess the degree of overlap between two or more sets of genes of interest using the R package euler which generates proportionally accurate sized venn diagrams unlike other common web-based venn diagram makers such as venny. For example, I’m curious to know the degree of overlap between upregulated DEGs in sgArid1a MC38 colon cancer cells and upregulated DEGs in sgArid1a B16F10 melanoma cell and will make a venn diagram below to get my answer.
#Filter for MC38 upregulated genes
MC38_upregulated <- MC38_edgeR_df %>%
filter(Log2fc >= .585, adj_p_value <= .05)
#Import in sgArid1a B16F10 RNA-seq data
B16F10_edgeR_df <- read.delim("sgArid1a_B16F10_edgeR.txt", header = TRUE, sep = "\t")
#replace single periods or double periods in column names with a single underscore
colnames(B16F10_edgeR_df) <- gsub("\\.|\\.\\.", "_", colnames(B16F10_edgeR_df))
#Rename column names in B16F10 edgeR dataframe
B16F10_edgeR_df <- rename(B16F10_edgeR_df, Gene = Annotation_Divergence, Log2fc = sgARID1A_vs_sgNT_Log2_Fold_Change, adj_p_value = sgARID1A_vs_sgNT_adj_p_value)
#Grab the gene name values preceding "|" edgeR from Gene column
B16F10_edgeR_df$Gene <- sub("\\|.*", "", B16F10_edgeR_df$Gene)
#Filter for upregulated B16 genes
B16_upregulated <- B16F10_edgeR_df %>%
filter(Log2fc >= .585, adj_p_value <= .05)
#Grab gene names in upregulated gene lists
MC38_upregulated <- MC38_upregulated[, "Gene"]
B16_upregulated <- B16_upregulated[, "Gene"]
#Find common and unique genes
common <- intersect(MC38_upregulated, B16_upregulated)
MC38_unique <- setdiff(MC38_upregulated, B16_upregulated)
B16_unique <- setdiff(B16_upregulated, MC38_upregulated)
# Count the number of common entries
num_common <- length(common)
#Eulerr ven diagram
library("eulerr")
fit <- euler(c("sgArid1a_MC38_up" = 184, "sgArid1a_B16F10_up" = 187,
"sgArid1a_MC38_up&sgArid1a_B16F10_up" = 65),
shape = "ellipse")
plot(fit, fills = c("skyblue", "darkgoldenrod1"), font = 9, quantities = TRUE)
From the results of the venn diagram, I can see that there are 65
commonly upregulated DEGs following ARID1A loss in the MC38 and B16F10
cancer cell lines. Since I’m interested in a phenotype of enhanced
anti-tumor immunity that’s observed in both sgArid1a tumor
models, these commonly upregulated genes can provided clues towards the
molecular mechanism underlying anti-tumor immunity in both models.
Session Info
For reproducibility, my session info is listed below outlining the version of R and package versions that I’ve using for this vignette.
print(sessionInfo())## R version 4.3.1 (2023-06-16 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 11 x64 (build 22631)
##
## Matrix products: default
##
##
## locale:
## [1] LC_COLLATE=English_United States.utf8
## [2] LC_CTYPE=English_United States.utf8
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.utf8
##
## time zone: America/Los_Angeles
## tzcode source: internal
##
## attached base packages:
## [1] grid stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] eulerr_7.0.0 clusterProfiler_4.8.3 msigdbr_7.5.1
## [4] RColorBrewer_1.1-3 viridis_0.6.4 viridisLite_0.4.2
## [7] ComplexHeatmap_2.16.0 ggrepel_0.9.4 scales_1.3.0
## [10] lubridate_1.9.3 forcats_1.0.0 stringr_1.5.1
## [13] dplyr_1.1.2 purrr_1.0.2 readr_2.1.4
## [16] tidyr_1.3.0 tibble_3.2.1 ggplot2_3.4.4
## [19] tidyverse_2.0.0 here_1.0.1 pacman_0.5.1
## [22] htmlwidgets_1.6.3 DiagrammeR_1.0.10
##
## loaded via a namespace (and not attached):
## [1] rstudioapi_0.15.0 jsonlite_1.8.8 shape_1.4.6
## [4] magrittr_2.0.3 magick_2.8.1 farver_2.1.1
## [7] rmarkdown_2.25 fs_1.6.3 GlobalOptions_0.1.2
## [10] zlibbioc_1.46.0 vctrs_0.6.3 memoise_2.0.1
## [13] RCurl_1.98-1.13 ggtree_3.8.2 htmltools_0.5.7
## [16] gridGraphics_0.5-1 sass_0.4.8 bslib_0.6.1
## [19] plyr_1.8.9 cachem_1.0.8 igraph_1.5.1
## [22] lifecycle_1.0.4 iterators_1.0.14 pkgconfig_2.0.3
## [25] gson_0.1.0 Matrix_1.6-4 R6_2.5.1
## [28] fastmap_1.1.1 GenomeInfoDbData_1.2.10 clue_0.3-65
## [31] digest_0.6.33 aplot_0.2.2 enrichplot_1.20.3
## [34] colorspace_2.1-0 patchwork_1.1.3 AnnotationDbi_1.62.2
## [37] S4Vectors_0.38.1 rprojroot_2.0.4 RSQLite_2.3.3
## [40] labeling_0.4.3 fansi_1.0.4 timechange_0.2.0
## [43] httr_1.4.7 polyclip_1.10-6 compiler_4.3.1
## [46] bit64_4.0.5 withr_2.5.2 doParallel_1.0.17
## [49] downloader_0.4 BiocParallel_1.34.2 DBI_1.1.3
## [52] highr_0.10 ggforce_0.4.1 MASS_7.3-60
## [55] rjson_0.2.21 HDO.db_0.99.1 tools_4.3.1
## [58] scatterpie_0.2.1 ape_5.7-1 glue_1.6.2
## [61] nlme_3.1-162 GOSemSim_2.26.1 polylabelr_0.2.0
## [64] shadowtext_0.1.2 cluster_2.1.4 reshape2_1.4.4
## [67] snow_0.4-4 fgsea_1.26.0 generics_0.1.3
## [70] gtable_0.3.4 tzdb_0.4.0 data.table_1.14.8
## [73] hms_1.1.3 tidygraph_1.2.3 utf8_1.2.3
## [76] XVector_0.40.0 BiocGenerics_0.46.0 foreach_1.5.2
## [79] pillar_1.9.0 yulab.utils_0.1.0 babelgene_22.9
## [82] circlize_0.4.15 splines_4.3.1 tweenr_2.0.2
## [85] treeio_1.24.3 lattice_0.21-8 bit_4.0.5
## [88] tidyselect_1.2.0 GO.db_3.17.0 Biostrings_2.68.1
## [91] knitr_1.45 gridExtra_2.3 IRanges_2.34.1
## [94] stats4_4.3.1 xfun_0.41 graphlayouts_1.0.2
## [97] Biobase_2.60.0 matrixStats_1.1.0 visNetwork_2.1.2
## [100] stringi_1.8.2 lazyeval_0.2.2 ggfun_0.1.3
## [103] yaml_2.3.7 evaluate_0.23 codetools_0.2-19
## [106] ggraph_2.1.0 qvalue_2.32.0 BiocManager_1.30.21
## [109] ggplotify_0.1.2 cli_3.6.1 munsell_0.5.0
## [112] jquerylib_0.1.4 Rcpp_1.0.11 GenomeInfoDb_1.36.4
## [115] png_0.1-8 parallel_4.3.1 ellipsis_0.3.2
## [118] blob_1.2.4 DOSE_3.26.2 bitops_1.0-7
## [121] tidytree_0.4.5 crayon_1.5.2 GetoptLong_1.0.5
## [124] rlang_1.1.1 cowplot_1.1.1 fastmatch_1.1-4
## [127] KEGGREST_1.40.1