Last updated: 2023-06-08 Code version: 503a830e86095d462d770554689962c078e8993b
library("tidyverse")
library("DESeq2")
library("clusterProfiler")This R markdown file showcases how I perform PCA and differential expression analysis on different comparison groups. Currently, it is only tested on Linux. The Galaxy workflow by Elisa Darkow was utilized to create featureCounts tables. They are available for download here. You can quickly unpack it in your working directory. Plus, the experimental design file is also downloadable from here.
read_table <- function(table_name, pat){
# Here I read table to dataframe and
# name the count column with sample ID.
# the input example:
# table_name = "featureCounts/SRR6898748.tabular"
# pat = ".tabular"
# extract sample ID from table file name
basename_ <- basename(table_name)
basename <- gsub(pat, "", basename_)
# add sample ID as column name
df <- read.table(table_name, sep = "\t", header = TRUE,
col.names = c("gene_id", basename),
stringsAsFactors = FALSE)
return(df)
}
get_counts <- function(store_dir, pat, count_filter=0){
# get a list of dataframes
tables <- list.files(store_dir, pattern=pat, full.names = TRUE)
dfs <- lapply(tables, read_table, pat)
# merge the list of dataframes to one datadframe
dfs_joined <- Reduce(function(...) full_join(..., by = "gene_id"), dfs)
# replace NA to 0
# NA are where the values are missing in one or more libraries
dfs_joined <- mutate_all(dfs_joined, funs(replace(., which(is.na(.)), 0)))
# here with count_filter, I can choose to filter
# the gene with low counts e.g. 10 reads
# this reduces the memory size by removing noisy low counts
# however, some novice biologist still wants the Genes
# with 0 counts for some reason.
dfs_filtered <- subset(dfs_joined, rowSums(subset(dfs_joined, select = 2:length(dfs_joined))) >= count_filter)
# replace rownames with column gene_id
dfs_filtered <- dfs_filtered %>% remove_rownames %>% column_to_rownames(var="gene_id")
counts <- as.matrix(dfs_filtered)
return(counts)
}
featureCounts_dir <- './featureCounts'
# merge counts of different samples into one dataframe
featureCounts <- get_counts(featureCounts_dir, '.csv')experimental_design <- './Experimental_design.csv'
# set all columns as factors
design <- read.table(experimental_design, sep = "\t", header = TRUE, row.names="sample_ID", colClasses='factor')
# change 'm' to 'male' and 'f' to 'female'
design$Sex <- str_replace(design$Sex, 'm', 'male') %>% str_replace(., 'f', 'female')
# order columns of count matrix according to rows of design data
featureCounts <- featureCounts[, design %>% rownames]
all(rownames(design) == colnames(featureCounts))## [1] TRUEmk_PCA <- function(dds_data, title, param_str, file_id, color="source", rev_PC1 = FALSE, rev_PC2 = FALSE, shape_title, shape_symbol){
# function to generate a PCA table
# input examples:
# dds_data = dds_RA_AFvsSR_TD
# title = 'RA, AF vs. SR, Thomas vs. Darkow'
# param_str = "Health.status..general."
# file_id = 'RA_AFvsSR_TD'
# transform raw counts with vst
vst_data <- vst(dds_data, blind=FALSE)
counts_vst <- vst_data %>% assay
# remove batch effects from vst transformed counts
counts_rb <- limma::removeBatchEffect(counts_vst, vst_data$source)
# save normalized counts after the batch removal
# var: file_id
# change rownames to colnames
counts_rb_ <- counts_rb %>% as.data.frame %>% rownames_to_column(var='ENSEMBL')
# run PCA
vst_rb <- vst_data
assay(vst_rb) <- counts_rb
# var: param_str
pcadata_vst_rb <- plotPCA(vst_rb, intgroup = c(param_str, color, 'Sex'),
returnData=TRUE)
# param_str %in% colnames(pcadata_vst_rb)
pcadata_vst_rb$group_gender <- paste(pcadata_vst_rb[,param_str], pcadata_vst_rb$Sex, sep=' ')
# save table of PCA parameters
# var: file_id
fn_PCA_table = paste("figures/pcadata_vst_rb_", file_id, ".csv", sep="")
write.csv(pcadata_vst_rb, file = fn_PCA_table, row.names=FALSE)
paste(fn_PCA_table, "is created.", sep=" ") %>% print
percentVar_vst <- round(100 * attr(pcadata_vst_rb, "percentVar"))
# var: param_str
if (rev_PC1){
PC1 = -pcadata_vst_rb$"PC1"
}else{PC1 = pcadata_vst_rb$"PC1"}
if (rev_PC2){
PC2 = -pcadata_vst_rb$"PC2"
}else{PC2 = pcadata_vst_rb$"PC2"}
ggplot(pcadata_vst_rb, aes_string(PC1, PC2, color=color, shape="group_gender")) +
geom_point(size=4) +
scale_shape_manual(values=shape_symbol)+
# geom_text(aes(label=name),hjust=0, vjust=0) +
xlab(paste0("PC1: ",percentVar_vst[1],"% variance")) +
ylab(paste0("PC2: ",percentVar_vst[2],"% variance")) +
# ggtitle(title) +
scale_colour_brewer(palette="Set1") +
labs(shape=shape_title, color = "Reference") +
theme(text = element_text(size=12))
# var: file_id
fn_PCA = paste("figures/PCA_vst_rb_", file_id, ".png", sep="")
ggsave(fn_PCA)
paste(fn_PCA, "is created.", sep=" ") %>% print
return(counts_rb_)
}
mk_DE <- function(dds, condition, factor_1, factor_2, file_id, counts_rb, padj_filter = 0.05){
# This function is to generate DE table
# input examples:
# dds = dds_deseq_LV_HFvsdornor_LSD
# condition = 'Health.status..general.'
# factor_1 = 'HF'
# factor_2 = 'Donor'
# file_id = 'LV_HFvsdornor_LSD'
# counts_rb = counts_rb_HFvsdornor_LSD
# calculate DE
# to account for noisy genes with low counts or high variability
contrast <- lfcShrink(dds=dds, contrast=c(condition, factor_1, factor_2), type='ashr')
contrast_lfc <- contrast[, c('baseMean', 'log2FoldChange', 'padj')] %>%
as.data.frame %>%
subset(., padj <= 0.05) %>% # filter out gene without significance
rownames_to_column(var="ENSEMBL")
contrast_rb_cnt <- merge(contrast_lfc, counts_rb, by = 'ENSEMBL')
counts_dds <- counts(dds) %>% as.data.frame %>% rownames_to_column(var='ENSEMBL')
contrast_dds_cnt <- merge(contrast_lfc, counts_dds, by = 'ENSEMBL')
# translate ids
translate_ids <- bitr(contrast_rb_cnt[, 'ENSEMBL'], fromType="ENSEMBL", toType=c("ENTREZID", "SYMBOL", "GENENAME"), OrgDb="org.Hs.eg.db", drop = FALSE)
contrast_rb_cnt_symbol <- merge(translate_ids, contrast_rb_cnt, by = 'ENSEMBL') %>%
.[order(.$padj),]
contrast_dds_cnt_symbol <- merge(translate_ids, contrast_dds_cnt, by = 'ENSEMBL') %>%
.[order(.$padj),]
fn_DE = paste("tables/DE_rbCount_", file_id, ".csv", sep="")
# var: file_id
write.csv(contrast_rb_cnt_symbol, file = fn_DE, row.names=FALSE)
print(paste(fn_DE, "is created.", sep=" "))
fn_DE_dds = paste("tables/DE_ddsCount_", file_id, ".csv", sep="")
# var: file_id
write.csv(contrast_dds_cnt_symbol, file = fn_DE_dds, row.names=FALSE)
print(paste(fn_DE_dds, "is created.", sep=" "))
}
# create directories for figures and tables
create_dir <- function(directory){
if (!dir.exists(directory)) {
# Create the directory
dir.create(directory)
print(paste("Directory", directory, "created."))
} else {
print(paste("Directory", directory, "already exists."))
}
}
create_dir('./figures')## [1] "Directory ./figures created."create_dir('./tables')## [1] "Directory ./tables created."# subset: LV, HF vs. donor, Liu vs. Schiano vs. Darkow
# Prepare dataframes of design and counts
design_LV_HFvsdornor_LSD <- subset(design, Sample.provenance.specific. %in% c('LV')
& Health.status..general. %in% c('Donor', 'HF')
& source %in% c('Schiano', 'Liu', 'Darkow')) %>%
droplevels
counts_raw_LV_HFvsdornor_LSD <- featureCounts[, rownames(design_LV_HFvsdornor_LSD)]
all(rownames(design_LV_HFvsdornor_LSD) == colnames(counts_raw_LV_HFvsdornor_LSD))## [1] TRUE# create DESeqDataSet From Matrix
dds_LV_HFvsdornor_LSD <- DESeqDataSetFromMatrix(countData = counts_raw_LV_HFvsdornor_LSD,
colData = design_LV_HFvsdornor_LSD,
design = ~ source + Health.status..general.)
file_id = "LV_HFvsdornor_LSD"
# make PCA
counts_rb_HFvsdornor_LSD <- mk_PCA(dds_data = dds_LV_HFvsdornor_LSD,
title = 'LV, HF vs. donor, Liu vs. Schiano vs. Darkow',
color = "source",
param_str = "Health.status..general.",
rev_PC2 = TRUE, rev_PC1 = TRUE,
file_id = file_id,
shape_title = "Health status and gender",
shape_symbol = c(15, 0, 17, 2))## [1] "figures/pcadata_vst_rb_LV_HFvsdornor_LSD.csv is created."## [1] "figures/PCA_vst_rb_LV_HFvsdornor_LSD.png is created."# make DE
dds_deseq_LV_HFvsdornor_LSD <- DESeq(dds_LV_HFvsdornor_LSD)
mk_DE(dds = dds_deseq_LV_HFvsdornor_LSD,
condition = 'Health.status..general.',
factor_1 = 'HF', factor_2 = 'Donor',
file_id = file_id, counts_rb = counts_rb_HFvsdornor_LSD)## [1] "tables/DE_rbCount_LV_HFvsdornor_LSD.csv is created."
## [1] "tables/DE_ddsCount_LV_HFvsdornor_LSD.csv is created."# subset: LV, DCM vs. donor, Liu vs. Schiano vs. Darkow
design_LV_DCMvsdonor_LSD <- subset(design,
Sample.provenance.specific. %in% c('LV')
& Health.status..specific. %in% c('Donor', 'DCM')
& source %in% c('Schiano', 'Liu', 'Darkow')) %>%
droplevels
# design_LV_DCMvsdonor_LSD %>% kable(digits=100)
counts_raw_LV_DCMvsdonor_LSD <- featureCounts[, rownames(design_LV_DCMvsdonor_LSD)]
all(rownames(design_LV_DCMvsdonor_LSD) == colnames(counts_raw_LV_DCMvsdonor_LSD))## [1] TRUE# create DESeqDataSet From Matrix
# var: Health.status..specific.
dds_LV_DCMvsdonor_LSD <- DESeqDataSetFromMatrix(countData = counts_raw_LV_DCMvsdonor_LSD,
colData = design_LV_DCMvsdonor_LSD,
design = ~ source + Health.status..specific.)
file_id = 'LV_DCMvsdonor_LSD'
# make PCA
counts_rb_LV_DCMvsdonor_LSD <- mk_PCA(dds_data = dds_LV_DCMvsdonor_LSD,
title = 'LV, DCM vs. donor, Liu vs. Schiano vs. Darkow',
param_str = "Health.status..specific.",
file_id = file_id, shape_title="Health status and gender",
rev_PC2 = FALSE, shape_symbol=c(19, 1, 15, 0))## [1] "figures/pcadata_vst_rb_LV_DCMvsdonor_LSD.csv is created."## [1] "figures/PCA_vst_rb_LV_DCMvsdonor_LSD.png is created."# make DE
dds_deseq_LV_DCMvsdonor_LSD <- DESeq(dds_LV_DCMvsdonor_LSD)
mk_DE(dds = dds_deseq_LV_DCMvsdonor_LSD,
condition = 'Health.status..specific.',
factor_1 = 'DCM', factor_2 = 'Donor',
file_id = file_id, counts_rb = counts_rb_LV_DCMvsdonor_LSD)## [1] "tables/DE_rbCount_LV_DCMvsdonor_LSD.csv is created."
## [1] "tables/DE_ddsCount_LV_DCMvsdonor_LSD.csv is created."# colnames(design)
# 'Sample.provenance.specific.''Sample.provenance..general.'
# 'Health.status..specific.''Health.status..general.''Sex''Age''source'
# subset: RA, AF vs. SR, Thomas vs. Darkow
design_RA_AFvsSR_TD <- subset(design,
Sample.provenance.specific. %in% c('RA')
& Health.status..general. %in% c('AF', 'SR')
& source %in% c('Thomas', 'Darkow')) %>%
droplevels
# design_RA_AFvsSR_TD %>% kable(digits=100)
counts_raw_RA_AFvsSR_TD <- featureCounts[, rownames(design_RA_AFvsSR_TD)]
all(rownames(design_RA_AFvsSR_TD) == colnames(counts_raw_RA_AFvsSR_TD))## [1] TRUE# create DESeqDataSet From Matrix
dds_RA_AFvsSR_TD <- DESeqDataSetFromMatrix(countData = counts_raw_RA_AFvsSR_TD,
colData = design_RA_AFvsSR_TD,
design = ~ source + Health.status..general.)
file_id = 'RA_AFvsSR_TD'
# make PCA
counts_rb_RA_AFvsSR_TD <- mk_PCA(dds_data = dds_RA_AFvsSR_TD,
title = 'RA, AF vs. SR, Thomas vs. Darkow',
param_str = "Health.status..general.",
file_id = file_id, shape_title="Health status and gender",
rev_PC2 = FALSE, shape_symbol=c(5, 10))## [1] "figures/pcadata_vst_rb_RA_AFvsSR_TD.csv is created."## [1] "figures/PCA_vst_rb_RA_AFvsSR_TD.png is created."# make DE
dds_deseq_RA_AFvsSR_TD <- DESeq(dds_RA_AFvsSR_TD)
mk_DE(dds = dds_deseq_RA_AFvsSR_TD,
condition = 'Health.status..general.',
factor_1 = 'AF', factor_2 = 'SR',
file_id = file_id, counts_rb = counts_rb_RA_AFvsSR_TD)## [1] "tables/DE_rbCount_RA_AFvsSR_TD.csv is created."
## [1] "tables/DE_ddsCount_RA_AFvsSR_TD.csv is created."# subset: RA, AF vs. CAD, Thomas vs. Darkow
design_RA_AFvsCAD_TD <- subset(design,
Sample.provenance.specific. %in% c('RA')
& Health.status..specific. %in% c('AF', 'CAD')
& source %in% c('Thomas', 'Darkow')) %>%
droplevels
counts_raw_RA_AFvsCAD_TD <- featureCounts[, rownames(design_RA_AFvsCAD_TD)]
all(rownames(design_RA_AFvsCAD_TD) == colnames(counts_raw_RA_AFvsCAD_TD))## [1] TRUE# create DESeqDataSet From Matrix
dds_RA_AFvsCAD_TD <- DESeqDataSetFromMatrix(countData = counts_raw_RA_AFvsCAD_TD,
colData = design_RA_AFvsCAD_TD,
design = ~ source + Health.status..specific.)
file_id = 'RA_AFvsCAD_TD'
# make PCA
counts_rb_RA_AFvsCAD_TD <- mk_PCA(dds_data = dds_RA_AFvsCAD_TD,
title = 'RA, AF vs. CAD, Thomas vs. Darkow',
param_str = "Health.status..specific.",
file_id = file_id, shape_title="Health status and gender",
rev_PC2 = FALSE, shape_symbol=c(5, 10))## [1] "figures/pcadata_vst_rb_RA_AFvsCAD_TD.csv is created."## [1] "figures/PCA_vst_rb_RA_AFvsCAD_TD.png is created."# make DE
dds_deseq_RA_AFvsCAD_TD <- DESeq(dds_RA_AFvsCAD_TD)
mk_DE(dds = dds_deseq_RA_AFvsCAD_TD,
condition = 'Health.status..specific.',
factor_1 = 'AF', factor_2 = 'CAD',
file_id = file_id, counts_rb = counts_rb_RA_AFvsCAD_TD)## [1] "tables/DE_rbCount_RA_AFvsCAD_TD.csv is created."
## [1] "tables/DE_ddsCount_RA_AFvsCAD_TD.csv is created."# subset: donor, atria vs. ventricles, Johnson vs. Darkow
design_donor_atriaVSventricles_JD <- subset(design,
Health.status..general. %in% c('Donor')
& Sample.provenance..general. %in% c('Atria', 'Ventricles')
& source %in% c('Johnson', 'Darkow')) %>%
droplevels
# design_donor_atriaVSventricles_JD %>% kable(digits=100)
counts_raw_donor_atriaVSventricles_JD <- featureCounts[, rownames(design_donor_atriaVSventricles_JD)]
all(rownames(design_donor_atriaVSventricles_JD) == colnames(counts_raw_donor_atriaVSventricles_JD))## [1] TRUE# create DESeqDataSet From Matrix
dds_donor_atriaVSventricles_JD <- DESeqDataSetFromMatrix(countData = counts_raw_donor_atriaVSventricles_JD,
colData = design_donor_atriaVSventricles_JD,
design = ~ source + Sample.provenance..general.)
file_id = 'donor_atriaVSventricles_JD'
# make PCA
counts_rb_donor_atriaVSventricles_JD <- mk_PCA(dds_data = dds_donor_atriaVSventricles_JD,
title = 'donor, atria vs. ventricles, Johnson vs. Darkow',
param_str = "Sample.provenance..general.",
file_id = file_id, rev_PC2 = TRUE,
shape_title = "Tissue provenance and gender",
shape_symbol = c(19, 1, 15, 0))## [1] "figures/pcadata_vst_rb_donor_atriaVSventricles_JD.csv is created."## [1] "figures/PCA_vst_rb_donor_atriaVSventricles_JD.png is created."# make DE
dds_deseq_donor_atriaVSventricles_JD <- DESeq(dds_donor_atriaVSventricles_JD)
mk_DE(dds = dds_deseq_donor_atriaVSventricles_JD,
condition = 'Sample.provenance..general.',
factor_1 = 'Atria', factor_2 = 'Ventricles',
file_id = file_id, counts_rb = counts_rb_donor_atriaVSventricles_JD,
padj_filter = 1) # Elisa wants keep all genes without filtering genes by counts## [1] "tables/DE_rbCount_donor_atriaVSventricles_JD.csv is created."
## [1] "tables/DE_ddsCount_donor_atriaVSventricles_JD.csv is created."# subset: donor, LA vs. LV, Johnson vs. Darkow
design_donor_LAvsLV_JD <- subset(design,
Health.status..general. %in% c('Donor')
& Sample.provenance.specific. %in% c('LA', 'LV')
& source %in% c('Johnson', 'Darkow')) %>%
droplevels
counts_raw_donor_LAvsLV_JD <- featureCounts[, rownames(design_donor_LAvsLV_JD)]
all(rownames(design_donor_LAvsLV_JD) == colnames(counts_raw_donor_LAvsLV_JD))## [1] TRUE# create DESeqDataSet From Matrix
dds_donor_LAvsLV_JD <- DESeqDataSetFromMatrix(countData = counts_raw_donor_LAvsLV_JD,
colData = design_donor_LAvsLV_JD,
design = ~ source + Sample.provenance.specific.)
file_id = 'donor_LAvsLV_JD'
# make PCA
counts_rb_donor_LAvsLV_JD <- mk_PCA(dds_data = dds_donor_LAvsLV_JD,
title = 'donor, LA vs. LV, Johnson vs. Darkow',
param_str = "Sample.provenance.specific.",
file_id = file_id, rev_PC2 = TRUE,
shape_title = "Sample provenance and gender",
shape_symbol = c(19, 1, 15, 0))## [1] "figures/pcadata_vst_rb_donor_LAvsLV_JD.csv is created."## [1] "figures/PCA_vst_rb_donor_LAvsLV_JD.png is created."# make DE
dds_deseq_donor_LAvsLV_JD <- DESeq(dds_donor_LAvsLV_JD)
mk_DE(dds = dds_deseq_donor_LAvsLV_JD,
condition = 'Sample.provenance.specific.',
factor_1 = 'LA', factor_2 = 'LV',
file_id = file_id, counts_rb = counts_rb_donor_LAvsLV_JD)## [1] "tables/DE_rbCount_donor_LAvsLV_JD.csv is created."
## [1] "tables/DE_ddsCount_donor_LAvsLV_JD.csv is created."# subset: donor, RA vs. LA, Johnson vs. Darkow
design_donor_RAvsLA_JD <- subset(design,
Health.status..general. %in% c('Donor')
& Sample.provenance.specific. %in% c('RA', 'LA')
& source %in% c('Johnson', 'Darkow')) %>%
droplevels
# design_donor_RAvsLA_JD %>% kable(digits=100)
counts_raw_donor_RAvsLA_JD <- featureCounts[, rownames(design_donor_RAvsLA_JD)]
all(rownames(design_donor_RAvsLA_JD) == colnames(counts_raw_donor_RAvsLA_JD))## [1] TRUE# create DESeqDataSet From Matrix
dds_donor_RAvsLA_JD <- DESeqDataSetFromMatrix(countData = counts_raw_donor_RAvsLA_JD,
colData = design_donor_RAvsLA_JD,
design = ~ source + Sample.provenance.specific.)
file_id = 'donor_RAvsLA_JD'
# make PCA
counts_rb_donor_RAvsLA_JD <- mk_PCA(dds_data = dds_donor_RAvsLA_JD,
title = 'donor, RA vs. LA, Johnson vs. Darkow',
param_str = "Sample.provenance.specific.",
file_id = file_id, rev_PC2 = TRUE,
shape_title = "Sample provenance and gender",
shape_symbol = c(19, 1, 15, 0))## [1] "figures/pcadata_vst_rb_donor_RAvsLA_JD.csv is created."## [1] "figures/PCA_vst_rb_donor_RAvsLA_JD.png is created."# make DE
dds_deseq_donor_RAvsLA_JD <- DESeq(dds_donor_RAvsLA_JD)
mk_DE(dds = dds_deseq_donor_RAvsLA_JD,
condition = 'Sample.provenance.specific.',
factor_1 = 'RA', factor_2 = 'LA',
file_id = file_id, counts_rb = counts_rb_donor_RAvsLA_JD)## [1] "tables/DE_rbCount_donor_RAvsLA_JD.csv is created."
## [1] "tables/DE_ddsCount_donor_RAvsLA_JD.csv is created."sessionInfo()## R version 4.2.2 (2022-10-31)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.6 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/libopenblasp-r0.2.20.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=de_DE.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=de_DE.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=de_DE.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=de_DE.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] clusterProfiler_4.4.4 DESeq2_1.36.0
## [3] SummarizedExperiment_1.26.1 Biobase_2.56.0
## [5] MatrixGenerics_1.8.1 matrixStats_0.62.0
## [7] GenomicRanges_1.48.0 GenomeInfoDb_1.32.4
## [9] IRanges_2.30.1 S4Vectors_0.34.0
## [11] BiocGenerics_0.42.0 forcats_0.5.2
## [13] stringr_1.4.1 dplyr_1.0.10
## [15] purrr_0.3.5 readr_2.1.3
## [17] tidyr_1.2.1 tibble_3.1.8
## [19] ggplot2_3.3.6 tidyverse_1.3.2
##
## loaded via a namespace (and not attached):
## [1] shadowtext_0.1.2 readxl_1.4.1 backports_1.4.1
## [4] fastmatch_1.1-3 systemfonts_1.0.4 plyr_1.8.7
## [7] igraph_1.3.5 lazyeval_0.2.2 splines_4.2.2
## [10] BiocParallel_1.30.4 digest_0.6.30 invgamma_1.1
## [13] yulab.utils_0.0.5 htmltools_0.5.3 GOSemSim_2.22.0
## [16] viridis_0.6.2 GO.db_3.15.0 SQUAREM_2021.1
## [19] fansi_1.0.3 magrittr_2.0.3 memoise_2.0.1
## [22] googlesheets4_1.0.1 limma_3.52.4 tzdb_0.3.0
## [25] Biostrings_2.64.1 annotate_1.74.0 graphlayouts_0.8.3
## [28] modelr_0.1.9 enrichplot_1.16.2 colorspace_2.0-3
## [31] blob_1.2.3 rvest_1.0.3 ggrepel_0.9.1
## [34] textshaping_0.3.6 haven_2.5.1 xfun_0.34
## [37] crayon_1.5.2 RCurl_1.98-1.9 jsonlite_1.8.2
## [40] scatterpie_0.1.8 genefilter_1.78.0 survival_3.4-0
## [43] ape_5.6-2 glue_1.6.2 polyclip_1.10-0
## [46] gtable_0.3.1 gargle_1.2.1 zlibbioc_1.42.0
## [49] XVector_0.36.0 DelayedArray_0.22.0 scales_1.2.1
## [52] DOSE_3.22.1 DBI_1.1.3 Rcpp_1.0.9
## [55] viridisLite_0.4.1 xtable_1.8-4 gridGraphics_0.5-1
## [58] tidytree_0.4.1 bit_4.0.4 truncnorm_1.0-8
## [61] httr_1.4.4 fgsea_1.22.0 RColorBrewer_1.1-3
## [64] ellipsis_0.3.2 pkgconfig_2.0.3 XML_3.99-0.11
## [67] farver_2.1.1 sass_0.4.2 dbplyr_2.2.1
## [70] locfit_1.5-9.6 utf8_1.2.2 labeling_0.4.2
## [73] ggplotify_0.1.0 tidyselect_1.2.0 rlang_1.0.6
## [76] reshape2_1.4.4 AnnotationDbi_1.58.0 munsell_0.5.0
## [79] cellranger_1.1.0 tools_4.2.2 cachem_1.0.6
## [82] downloader_0.4 cli_3.4.1 generics_0.1.3
## [85] RSQLite_2.2.18 broom_1.0.1 evaluate_0.17
## [88] fastmap_1.1.0 ragg_1.2.5 yaml_2.3.6
## [91] ggtree_3.4.4 org.Hs.eg.db_3.15.0 knitr_1.40
## [94] bit64_4.0.5 fs_1.5.2 tidygraph_1.2.2
## [97] KEGGREST_1.36.3 ggraph_2.1.0 nlme_3.1-160
## [100] aplot_0.1.8 DO.db_2.9 xml2_1.3.3
## [103] compiler_4.2.2 png_0.1-7 treeio_1.20.2
## [106] reprex_2.0.2 tweenr_2.0.2 geneplotter_1.74.0
## [109] bslib_0.4.0 stringi_1.7.8 lattice_0.20-45
## [112] Matrix_1.5-1 vctrs_0.4.2 pillar_1.8.1
## [115] lifecycle_1.0.3 jquerylib_0.1.4 irlba_2.3.5.1
## [118] data.table_1.14.4 bitops_1.0-7 patchwork_1.1.2
## [121] qvalue_2.28.0 R6_2.5.1 gridExtra_2.3
## [124] codetools_0.2-18 MASS_7.3-58.1 assertthat_0.2.1
## [127] withr_2.5.0 GenomeInfoDbData_1.2.8 parallel_4.2.2
## [130] hms_1.1.2 grid_4.2.2 ggfun_0.0.7
## [133] rmarkdown_2.17 ashr_2.2-54 googledrive_2.0.0
## [136] mixsqp_0.3-43 ggforce_0.4.1 lubridate_1.8.0