R/Bioconductorでマイクロアレイ解析
大学で今日やったハンズオン講習会は学内環境の問題で上手く行かなかった。勿体無いので資料を公開します。 もっとこうした方が良いとか、ここが間違っているとかあったら是非コメントしてください。 初投稿なので勝手が分からない...
Microarray analysis ~Is microarray analysis outdated?~
Let’s analyse microarray with R/Bioconductor. Compared with RNA-seq, microarray requires less CPU and memory consumption. We no longer need to pay for expensive softwares such as GeneSpring to get gene expression data and conduct differential gene expression analysis. R/Bioconductor is an open-source and freely available repository for bioinformatics. You can almost do anything you want with R/Bioconductor.
I used the workflow below as reference http://bioconductor.org/packages/release/workflows/html/maEndToEnd.html.
FYI, the latest version of Bioconductor 3.10 is now availablehttps://www.bioconductor.org.
# 20191102 handson trainning for microarray analysis at YCU, Fuku-ura
# Affymetrics microarray
# download and read CEL files, RMA normalisation, quality check (PCA, clustering and so on), DEA, visualisation such as MA plot, volcano plot and heatmap.
# about the data
# WT and Bcl6 cKO (probably Itgax-cre) mice splenic DCs(cDC1, cDC2, pDC)
# Kenneth Murphy's lab, unpublished data
# we are going to compare WT vs Bcl6cKO cDC1 this time
# accession number is GSE135904 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE135904
# platform [MoGene-1_0-st] Affymetrix Mouse Gene 1.0 ST Array [transcript (gene) version]
##############
# install and load required packages
library(BiocManager)
# BiocManager::install(c("mogene10sttranscriptcluster.db", "pd.mogene.1.0.st.v1", GEOquery"))
library(GEOquery)
library(mogene10sttranscriptcluster.db)
library(pd.mogene.1.0.st.v1)
library(oligo)
library(limma)
# general packeges
# install.packages()
library(tidyverse)
library(gplots)
library(ggrepel)
library(viridis)
library(RColorBrewer)
############### # download the raw data from GEO # download manually from https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE135904 # or # automatically download # using GEOquery package # getGEOSuppFiles('GSE135904', fetch_files = TRUE) # download list.files("GSE135904/GSE135904_RAW/")
## [1] "GSM4037643_WT_cDC1_rep1.CEL.gz" "GSM4037644_WT_cDC1_rep2.CEL.gz"
## [3] "GSM4037645_WT_cDC1_rep3.CEL.gz" "GSM4037646_cKO_cDC1_rep1.CEL.gz"
## [5] "GSM4037647_cKO_cDC1_rep2.CEL.gz" "GSM4037648_cKO_cDC1_rep3.CEL.gz"
## [7] "GSM4037649_cKO_cDC1_rep4.CEL.gz" "GSM4037650_WT_cDC2_rep1.CEL.gz"
## [9] "GSM4037651_WT_cDC2_rep2.CEL.gz" "GSM4037652_WT_cDC2_rep3.CEL.gz"
## [11] "GSM4037653_cKO_cDC2_rep1.CEL.gz" "GSM4037654_cKO_cDC2_rep2.CEL.gz"
## [13] "GSM4037655_cKO_cDC2_rep3.CEL.gz" "GSM4037656_cKO_cDC2_rep4.CEL.gz"
## [15] "GSM4037657_cKO_cDC2_rep5.CEL.gz" "GSM4037658_WT_pDC_rep1.CEL.gz"
## [17] "GSM4037659_WT_pDC_rep2.CEL.gz" "GSM4037660_WT_pDC_rep3.CEL.gz"
## [19] "GSM4037661_cKO_pDC_rep1.CEL.gz" "GSM4037662_cKO_pDC_rep2.CEL.gz"
## [21] "GSM4037663_cKO_pDC_rep3.CEL.gz" "GSM4037664_cKO_pDC_rep4.CEL.gz"
## [23] "GSM4037665_cKO_pDC_rep5.CEL.gz"
# read CEL files # extract only cDC1 data celFiles <- list.celfiles("GSE135904/GSE135904_RAW/", listGzipped = T, full.names = TRUE) celFiles2 <- celFiles %>% str_subset(pattern = "cDC1") # extract which contains "cDC1" in its name affyRaw <- oligo::read.celfiles(celFiles2) # reading CEL files
## Platform design info loaded.
## Reading in : GSE135904/GSE135904_RAW//GSM4037643_WT_cDC1_rep1.CEL.gz
## Reading in : GSE135904/GSE135904_RAW//GSM4037644_WT_cDC1_rep2.CEL.gz
## Reading in : GSE135904/GSE135904_RAW//GSM4037645_WT_cDC1_rep3.CEL.gz
## Reading in : GSE135904/GSE135904_RAW//GSM4037646_cKO_cDC1_rep1.CEL.gz
## Reading in : GSE135904/GSE135904_RAW//GSM4037647_cKO_cDC1_rep2.CEL.gz
## Reading in : GSE135904/GSE135904_RAW//GSM4037648_cKO_cDC1_rep3.CEL.gz
## Reading in : GSE135904/GSE135904_RAW//GSM4037649_cKO_cDC1_rep4.CEL.gz
# next step is normalisation # RMA - Robust Multichip Average algorithm # Robust Multichip Average preprocessing methodology. # This strategy allows background subtraction, quantile normalization and summarization (via median-polish). # RMA method is now the most popular normalisation methods for Affymetrics microarray. There are other methods like MAS5 and gcrma. eset <- oligo::rma(affyRaw) # RMA
## Background correcting
## Normalizing
## Calculating Expression
#################### # quality check # boxplot par(mfrow = c(1,2)) # separate plot panel oligo::boxplot(affyRaw, target = "core", las = 3, main = "raw data") # before RMA oligo::boxplot(eset, las = 3, main = "normalised data (RMA method)") # after RMA
# density plot oligo::hist(affyRaw, target = "core", main = "raw data")
oligo::hist(eset, main = "normalised data (RMA method)")
# output an expression matrix Biobase::write.exprs(eset, file = "cDC1_gene_expression.txt") # change column names for simplicity cDC1_gene_expression <- Biobase::exprs(eset) colnames(cDC1_gene_expression) <- c("WT_1", "WT_2", "WT_3", "cKO_1", "cKO_2", "cKO_3", "cKO_4") # perform hierarchical clustering # several methods dist <- as.dist(1 - cor(cDC1_gene_expression, method = "pearson")) # calculate Pearson's distance cluster_ward <- hclust(dist, method = "ward.D2") # ward method cluster_ave <- hclust(dist, method = "average") # group average method # par(mfrow = c(1,2)) plot(cluster_ward, main = "ward", hang = -1)
plot(cluster_ave, main = "average", hang = -1) # any differences b/w these?
# PCA principal component analysis # R uses `prcomp` function for PCA # generally speaking, an expression matrix has genes in its rows and samples in columns # so we need to transpose rows and columns # by using t() function in R PCA <- prcomp(t(cDC1_gene_expression), scale = FALSE) # data already normalised so that "scale = FALSE" percentVar <- round(100*PCA$sdev^2/sum(PCA$sdev^2),1) # calculate stats, these codes copied from somewhere on the Internet sd_ratio <- sqrt(percentVar[2] / percentVar[1]) is(PCA)
## [1] "prcomp"
is(PCA$x[, 1:2])
## [1] "matrix" "array" "mMatrix"
## [4] "ff_or_matrix" "structure" "vector"
## [7] "vector_OR_factor"
PCA.df <- PCA$x[, 1:2] # get PC1 and PC2 PCA.df
## PC1 PC2
## WT_1 -36.674123 -41.688626
## WT_2 49.474829 -29.495301
## WT_3 -5.812515 4.905374
## cKO_1 -19.223942 -7.888048
## cKO_2 35.817169 13.070479
## cKO_3 -19.262162 26.568017
## cKO_4 -4.319256 34.528105
# add phenodata column PCA.df <- transform(PCA.df, Pheno = c("WT", "WT", "WT", "cKO", "cKO", "cKO", "cKO")) PCA.df
## PC1 PC2 Pheno
## WT_1 -36.674123 -41.688626 WT
## WT_2 49.474829 -29.495301 WT
## WT_3 -5.812515 4.905374 WT
## cKO_1 -19.223942 -7.888048 cKO
## cKO_2 35.817169 13.070479 cKO
## cKO_3 -19.262162 26.568017 cKO
## cKO_4 -4.319256 34.528105 cKO
# plot using ggplot2 (better than base plot) ggplot(as.data.frame(PCA.df), aes(PC1, PC2)) + xlab(paste0("PC1, VarExp: ", percentVar[1], "%")) + ylab(paste0("PC2, VarExp: ", percentVar[2], "%")) + geom_point() + theme_bw() + geom_text_repel(label = colnames(cDC1_gene_expression))+ coord_fixed(ratio = 1)
# by colour ggplot(as.data.frame(PCA.df), aes(PC1, PC2, colour = Pheno)) + xlab(paste0("PC1, VarExp: ", percentVar[1], "%")) + ylab(paste0("PC2, VarExp: ", percentVar[2], "%")) + geom_point() + theme_bw() + geom_text_repel(label = colnames(cDC1_gene_expression)) + coord_fixed(ratio = 1)
############# # probe annotation # we need to know gene name instead of probeID # just one example here mogene10sttranscriptcluster.db
## ChipDb object:
## | DBSCHEMAVERSION: 2.1
## | Db type: ChipDb
## | Supporting package: AnnotationDbi
## | DBSCHEMA: MOUSECHIP_DB
## | ORGANISM: Mus musculus
## | SPECIES: Mouse
## | MANUFACTURER: Affymetrix
## | CHIPNAME: mogene10
## | MANUFACTURERURL: http://www.affymetrix.com
## | EGSOURCEDATE: 2017-Mar29
## | EGSOURCENAME: Entrez Gene
## | EGSOURCEURL: ftp://ftp.ncbi.nlm.nih.gov/gene/DATA
## | CENTRALID: ENTREZID
## | TAXID: 10090
## | GOSOURCENAME: Gene Ontology
## | GOSOURCEURL: ftp://ftp.geneontology.org/pub/go/godatabase/archive/latest-lite/
## | GOSOURCEDATE: 2017-Mar29
## | GOEGSOURCEDATE: 2017-Mar29
## | GOEGSOURCENAME: Entrez Gene
## | GOEGSOURCEURL: ftp://ftp.ncbi.nlm.nih.gov/gene/DATA
## | KEGGSOURCENAME: KEGG GENOME
## | KEGGSOURCEURL: ftp://ftp.genome.jp/pub/kegg/genomes
## | KEGGSOURCEDATE: 2011-Mar15
## | GPSOURCENAME: UCSC Genome Bioinformatics (Mus musculus)
## | GPSOURCEURL:
## | GPSOURCEDATE: 2017-Aug8
## | ENSOURCEDATE: 2017-Mar29
## | ENSOURCENAME: Ensembl
## | ENSOURCEURL: ftp://ftp.ensembl.org/pub/current_fasta
## | UPSOURCENAME: Uniprot
## | UPSOURCEURL: http://www.UniProt.org/
## | UPSOURCEDATE: Tue Sep 26 18:46:27 2017
##
## Please see: help('select') for usage information
pd.mogene.1.0.st.v1
## Class........: AffyGenePDInfo
## Manufacturer.: Affymetrix
## Genome Build.: MM10
## Chip Geometry: 1050 rows x 1050 columns
AnnotationDbi::columns(mogene10sttranscriptcluster.db)
## [1] "ACCNUM" "ALIAS" "ENSEMBL" "ENSEMBLPROT"
## [5] "ENSEMBLTRANS" "ENTREZID" "ENZYME" "EVIDENCE"
## [9] "EVIDENCEALL" "GENENAME" "GO" "GOALL"
## [13] "IPI" "MGI" "ONTOLOGY" "ONTOLOGYALL"
## [17] "PATH" "PFAM" "PMID" "PROBEID"
## [21] "PROSITE" "REFSEQ" "SYMBOL" "UNIGENE"
## [25] "UNIPROT"
key <- AnnotationDbi::keys(mogene10sttranscriptcluster.db, keytype = "PROBEID") annot <- AnnotationDbi::select(mogene10sttranscriptcluster.db, keys = key, columns = c("SYMBOL", "GENENAME"), keytype = "PROBEID")
## 'select()' returned 1:many mapping between keys and columns
# you can see many probes not assigned to any genes # we remove these probes from further analysis here annot_NA <- dplyr::filter(annot, is.na(SYMBOL)) annot <- dplyr::filter(annot, !is.na(SYMBOL)) # you may also see some probes assigned to multiple genes # remove multiple assigned probes too anno_grouped <- group_by(annot, PROBEID) anno_summarized <- dplyr::summarize(anno_grouped, no_of_matches = n_distinct(SYMBOL)) tail(anno_summarized)
## # A tibble: 6 x 2
## PROBEID no_of_matches
## <chr> <int>
## 1 10608608 37
## 2 10608613 21
## 3 10608615 23
## 4 10608625 12
## 5 10608628 21
## 6 10608630 35
anno_multimap <- filter(anno_summarized, no_of_matches > 1) head(anno_multimap)
## # A tibble: 6 x 2
## PROBEID no_of_matches
## <chr> <int>
## 1 10344713 2
## 2 10344719 2
## 3 10344741 3
## 4 10345030 2
## 5 10345079 3
## 6 10345083 3
# convert matrix to tibble # rowname (PROBEID) to 1st column cDC1_gene_expression_tb <- as.data.frame(cDC1_gene_expression) %>% tibble::rownames_to_column(var = "PROBEID") # remove unmapped or multi-mapped probes from expression data # %in% is an operater meaning if an element belongs to a vector cDC1_gene_expression_filtered <- cDC1_gene_expression_tb %>% dplyr::filter(!(PROBEID %in% c(anno_multimap$PROBEID, annot_NA$PROBEID))) cDC1_gene_expression_filtered <- cDC1_gene_expression_tb %>% dplyr::filter(!(PROBEID %in% c(anno_multimap$PROBEID, annot_NA$PROBEID))) # we conduct DEA without unmapped or multi-mapped probes str(cDC1_gene_expression)
## num [1:35556, 1:7] 12.51 5.19 10.86 9.36 2.54 ...
## - attr(*, "dimnames")=List of 2
## ..$ : chr [1:35556] "10338001" "10338002" "10338003" "10338004" ...
## ..$ : chr [1:7] "WT_1" "WT_2" "WT_3" "cKO_1" ...
str(cDC1_gene_expression_filtered)
## 'data.frame': 22098 obs. of 8 variables:
## $ PROBEID: chr "10344624" "10344633" "10344637" "10344653" ...
## $ WT_1 : num 7.72 8.15 7.25 3.34 6.77 ...
## $ WT_2 : num 7.03 8.65 7.86 3.45 6.98 ...
## $ WT_3 : num 7.72 8.63 7.38 3.34 6.85 ...
## $ cKO_1 : num 7.52 8.24 7.34 3.34 7.27 ...
## $ cKO_2 : num 7.83 8.4 7.64 3.4 6.87 ...
## $ cKO_3 : num 8.01 8.86 7.34 3.31 6.98 ...
## $ cKO_4 : num 7.8 8.7 7.28 3.14 7.02 ...
# convert data type for limma cDC1_gene_expression <- cDC1_gene_expression_filtered %>% tibble::column_to_rownames(var = "PROBEID") str(cDC1_gene_expression)
## 'data.frame': 22098 obs. of 7 variables:
## $ WT_1 : num 7.72 8.15 7.25 3.34 6.77 ...
## $ WT_2 : num 7.03 8.65 7.86 3.45 6.98 ...
## $ WT_3 : num 7.72 8.63 7.38 3.34 6.85 ...
## $ cKO_1: num 7.52 8.24 7.34 3.34 7.27 ...
## $ cKO_2: num 7.83 8.4 7.64 3.4 6.87 ...
## $ cKO_3: num 8.01 8.86 7.34 3.31 6.98 ...
## $ cKO_4: num 7.8 8.7 7.28 3.14 7.02 ...
################# # differential gene expression analysis by limma package # make contrast.matrix groups <- factor(c(0, 0, 0, 1, 1, 1, 1)) groups
## [1] 0 0 0 1 1 1 1
## Levels: 0 1
design <- model.matrix(~ 0 + groups) design
## groups0 groups1
## 1 1 0
## 2 1 0
## 3 1 0
## 4 0 1
## 5 0 1
## 6 0 1
## 7 0 1
## attr(,"assign")
## [1] 1 1
## attr(,"contrasts")
## attr(,"contrasts")$groups
## [1] "contr.treatment"
colnames(design) <- c("WT", "cKO") design
## WT cKO
## 1 1 0
## 2 1 0
## 3 1 0
## 4 0 1
## 5 0 1
## 6 0 1
## 7 0 1
## attr(,"assign")
## [1] 1 1
## attr(,"contrasts")
## attr(,"contrasts")$groups
## [1] "contr.treatment"
contrast.matrix <- makeContrasts(WTvscKO = WT - cKO, levels = design) contrast.matrix
## Contrasts
## Levels WTvscKO
## WT 1
## cKO -1
# fit to linear model fit <- lmFit(cDC1_gene_expression, design) fit2 <- contrasts.fit(fit, contrast.matrix) fit2 <- eBayes(fit2, proportion = 0.01) topTable(fit2)
## logFC AveExpr t P.Value adj.P.Val B
## 10534493 -6.550408 8.840045 -87.20756 3.624331e-11 8.009047e-07 8.136083
## 10364038 -3.300669 7.288904 -42.67066 3.642851e-09 4.024986e-05 7.658134
## 10576807 -4.827990 7.448193 -33.10386 1.865703e-08 1.374276e-04 7.288831
## 10498367 -2.802876 6.233213 -30.18603 3.374964e-08 1.551310e-04 7.114609
## 10493995 -4.152279 11.000441 -30.00204 3.510070e-08 1.551310e-04 7.102219
## 10548314 -3.328394 8.468356 -26.82584 7.195771e-08 2.650203e-04 6.855729
## 10539135 -4.969217 8.079582 -23.59303 1.637443e-07 5.169173e-04 6.524127
## 10576784 -2.623773 7.948250 -22.97149 1.942206e-07 5.364859e-04 6.448320
## 10362803 -2.181507 8.483018 -22.33110 2.326781e-07 5.713024e-04 6.365395
## 10573939 -2.403545 8.702084 -17.91644 9.463242e-07 1.880787e-03 5.624748
top250 <- topTable(fit2, adjust="fdr", sort.by="B", number=250) # summmary and visualisation summary(decideTests(fit2))
## WTvscKO
## Down 63
## NotSig 22032
## Up 3
limma::plotMA(fit2)
limma::plotSA(fit2)
limma::volcanoplot(fit2)
# save limma's result # convert again # There must be better way # PROBEID 1st column # tibbleのrownames_to_columnという便利な関数がある limma_result <- topTable(fit2, sort.by = "logFC", adjust.method = "fdr", number = Inf) %>% tibble::rownames_to_column(var = "PROBEID") # time to join # use left_join to combine limma result and probe annotation data. limma_result_gene <- dplyr::left_join(limma_result, annot, by = "PROBEID") write.table(limma_result_gene, file = "limma_result_gene.txt", sep = "\t", row.names = F) # multiple test # adjusted P value < 0.05 # as differential expressed genes limma_DEG <- limma_result_gene %>% dplyr::filter(adj.P.Val < 0.05) # finish
######## # visualisation # add gene symbol cDC1_gene_expression_filtered_gene <- cDC1_gene_expression_filtered %>% dplyr::left_join(annot, by = "PROBEID") Zeb2 <- filter(cDC1_gene_expression_filtered_gene, SYMBOL == "Zeb2") %>% select(2:8) barplot(as.matrix(Zeb2), xlab = "sample", ylab = "log2 expression", main = "Zeb2", las = 3)
Ccl24 <- filter(cDC1_gene_expression_filtered_gene, SYMBOL == "Ccl24") %>% select(2:8) barplot(as.matrix(Ccl24), xlab = "sample", ylab = "log2 expression", main = "Ccl24", las = 3)
# limma_result_merge <- merge(limma_result, anno_summarized, by = "PROBEID", all.x = F) # let's use ggplot2 # MA plot ma <- ggplot(data = limma_result, aes(x = AveExpr, y = logFC)) + geom_point(size = 0.4) + geom_point(data = limma_DEG, colour = "red", size = 0.7) + theme_classic() ma
# rabel gene symbol by ggrepel ma + geom_text_repel(data = limma_DEG, aes(label = SYMBOL))
# volcano plot vp <- ggplot(data = limma_result, aes(x = logFC, y = -log10(P.Value))) + geom_point(size = 0.3) + geom_point(data = limma_DEG, colour = "red", size = 0.3) + theme_classic() vp
vp + geom_text_repel(data = limma_DEG, aes(label = SYMBOL))
# heatmap # gplots::heatmap.2() heatmap_DEG <- cDC1_gene_expression_filtered_gene %>% filter(SYMBOL %in% limma_DEG$SYMBOL) %>% select(SYMBOL, everything()) %>% select(1,3:9) %>% distinct(SYMBOL, .keep_all = TRUE) %>% tibble::column_to_rownames(var = "SYMBOL") heatmap.2(as.matrix(heatmap_DEG), col = brewer.pal(11,"RdBu"), trace = "none", cexCol = 0.9, cexRow = 0.5, scale = "row")
heatmap.2(as.matrix(heatmap_DEG), col = brewer.pal(9,"YlGnBu"), trace = "none", cexCol = 0.9, cexRow = 0.5, scale = "row")
heatmap.2(as.matrix(heatmap_DEG), col = bluered(256), trace = "none", cexCol = 0.9, cexRow = 0.5, scale = "row")
Seurat <- colorRampPalette(c("#FF00FF", "#000000", "#FFFF00")) heatmap.2(as.matrix(heatmap_DEG), col = Seurat, trace = "none", cexCol = 0.9, cexRow = 0.5, scale = "row")
sessionInfo()
## R version 3.6.0 (2019-04-26)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Mojave 10.14.4
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
##
## locale:
## [1] ja_JP.UTF-8/ja_JP.UTF-8/ja_JP.UTF-8/C/ja_JP.UTF-8/ja_JP.UTF-8
##
## attached base packages:
## [1] stats4 parallel stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] RColorBrewer_1.1-2
## [2] viridis_0.5.1
## [3] viridisLite_0.3.0
## [4] ggrepel_0.8.1
## [5] gplots_3.0.1.1
## [6] forcats_0.4.0
## [7] stringr_1.4.0
## [8] dplyr_0.8.3
## [9] purrr_0.3.3
## [10] readr_1.3.1
## [11] tidyr_1.0.0
## [12] tibble_2.1.3
## [13] ggplot2_3.2.1
## [14] tidyverse_1.2.1
## [15] limma_3.40.6
## [16] pd.mogene.1.0.st.v1_3.14.1
## [17] DBI_1.0.0
## [18] oligo_1.48.0
## [19] oligoClasses_1.46.0
## [20] RSQLite_2.1.2
## [21] Biostrings_2.52.0
## [22] XVector_0.24.0
## [23] mogene10sttranscriptcluster.db_8.7.0
## [24] org.Mm.eg.db_3.8.2
## [25] AnnotationDbi_1.46.1
## [26] IRanges_2.18.3
## [27] S4Vectors_0.22.1
## [28] GEOquery_2.52.0
## [29] Biobase_2.44.0
## [30] BiocGenerics_0.30.0
## [31] BiocManager_1.30.9
##
## loaded via a namespace (and not attached):
## [1] nlme_3.1-141 bitops_1.0-6
## [3] matrixStats_0.55.0 lubridate_1.7.4
## [5] bit64_0.9-7 httr_1.4.1
## [7] GenomeInfoDb_1.20.0 tools_3.6.0
## [9] backports_1.1.5 utf8_1.1.4
## [11] R6_2.4.0 affyio_1.54.0
## [13] KernSmooth_2.23-16 lazyeval_0.2.2
## [15] colorspace_1.4-1 withr_2.1.2
## [17] gridExtra_2.3 tidyselect_0.2.5
## [19] bit_1.1-14 compiler_3.6.0
## [21] preprocessCore_1.46.0 cli_1.1.0
## [23] rvest_0.3.4 xml2_1.2.2
## [25] DelayedArray_0.10.0 labeling_0.3
## [27] caTools_1.17.1.2 scales_1.0.0
## [29] digest_0.6.22 rmarkdown_1.16
## [31] pkgconfig_2.0.3 htmltools_0.4.0
## [33] readxl_1.3.1 rlang_0.4.1
## [35] rstudioapi_0.10 generics_0.0.2
## [37] jsonlite_1.6 gtools_3.8.1
## [39] BiocParallel_1.18.1 RCurl_1.95-4.12
## [41] magrittr_1.5 GenomeInfoDbData_1.2.1
## [43] Matrix_1.2-17 fansi_0.4.0
## [45] Rcpp_1.0.2 munsell_0.5.0
## [47] lifecycle_0.1.0 stringi_1.4.3
## [49] yaml_2.2.0 SummarizedExperiment_1.14.1
## [51] zlibbioc_1.30.0 grid_3.6.0
## [53] affxparser_1.56.0 blob_1.2.0
## [55] gdata_2.18.0 crayon_1.3.4
## [57] lattice_0.20-38 haven_2.1.1
## [59] splines_3.6.0 hms_0.5.2
## [61] zeallot_0.1.0 knitr_1.25
## [63] pillar_1.4.2 GenomicRanges_1.36.1
## [65] codetools_0.2-16 glue_1.3.1
## [67] evaluate_0.14 modelr_0.1.5
## [69] vctrs_0.2.0 foreach_1.4.7
## [71] cellranger_1.1.0 gtable_0.3.0
## [73] assertthat_0.2.1 xfun_0.10
## [75] broom_0.5.2 ff_2.2-14
## [77] iterators_1.0.12 memoise_1.1.0
