Characterising newly made gene sets
The DEE2 Gene Signatures Group
2020-10-19
Source: https://github.com/markziemann/dee2_gene_signatures
Background
In this analysis we will be characterising the new gene sets.
Libraries
suppressPackageStartupMessages({
library("getDEE2")
library("mitch")
library("triwise")
library("dplyr")
library("gplots")
library("reshape2")
library("network")
library("eulerr")
})Reference gene list
This is the background - all genes in the "universe" according to Ensembl version 90, which is the same annotation set used by DEE2. You can see that there are a large number of non-protein coding genes.
# get universe of gene names and biotypes
if (! file.exists("hs.gtf.gz") ) {
download.file("ftp://ftp.ensembl.org/pub/release-90/gtf/homo_sapiens/Homo_sapiens.GRCh38.90.gtf.gz",destfile= "hs.gtf.gz")
}
g<-read.table("hs.gtf.gz",sep="\t")
g <- g[grep("gene",g$V3),9]
universe <- sapply(strsplit(g," "),"[[",6)
universe <- gsub(";","",universe)
biotypes <- sapply(strsplit(g," "),"[[",10)
biotypes <- gsub("_"," ",biotypes)
biotypes <- gsub(";","",biotypes)
biotypes_df <- data.frame(universe,biotypes)
universe <- unique(universe)
mytable <- table(biotypes)
mytable <- mytable[order(mytable)]
mytable <- mytable[which(mytable>100)]
par(mar=c(3,13,1,1)); barplot(mytable,horiz=TRUE,las=2,cex.names = 0.6,cex.axis = 0.6,
main="number of genes in each biotype class",xlim = c(0,20000)) ;grid()
mytable <- mytable/sum(mytable)*100
par(mar=c(3,13,1,1)); barplot(mytable,horiz=TRUE,las=2,cex.names = 0.6,cex.axis = 0.6,
main="proportion (%) of genes in each biotype class" ,xlim=c(0,40)) ;grid()
# Take note how many transcripts are protein coding versus non-protein coding
paste("protein coding:",length(which(biotypes=="protein coding")))## [1] "protein coding: 19847"paste("non-protein coding:",length(which(biotypes!="protein coding")))## [1] "non-protein coding: 38455"Curated gene sets
Reactome
Here I download the current Reactome gene set library. I would like to see how many Ensembl genes have some sort of annotated function. Breaking it down into protein coding and non-protein coding, we can see that non-protein coding genes are severely underrepresented in Reactome. Reactome has 2400 sets and a total of 11193 genes as of 29/Sep/2020.
#o
download.file("https://reactome.org/download/current/ReactomePathways.gmt.zip", destfile="ReactomePathways.gmt.zip")
#unzip("ReactomePathways.gmt.zip")
reactome <- gmt_import("ReactomePathways.gmt")
reactome_genes <- unique(unname(unlist(reactome)))
length(reactome)## [1] 2408length(reactome_genes)## [1] 11194v1 <- list("Reactome"=reactome_genes, "Ensembl universe"=universe)
library("eulerr")
plot(euler(v1),quantities = TRUE)
length(reactome_genes)/length(universe)*100## [1] 19.76063prot <- biotypes_df[which(biotypes_df$biotypes=="protein coding"),1]
prot <- unique(prot)
nprot <- biotypes_df[which(biotypes_df$biotypes!="protein coding"),1]
nprot <- unique(nprot)
intersect(prot,nprot)## [1] "RGS5" "TMEM247" "CYB561D2" "CDR1" "KBTBD11-OT1"
## [6] "MAL2" "SPATA13" "SFTA3" "GOLGA8M" "ELFN2"nprot <- setdiff(nprot,prot)
v1 <- list("Reactome"=reactome_genes, "protein coding"=prot, "non-protein coding"=nprot)
plot(euler(v1),quantities = TRUE)
Gene Ontology
Now let's do the same with GO sets downloaded from Ensembl biomart. GO has 12264 sets and 8412 genes as of 26/Sep/2020.
go <- read.table("biomart_2020-09-26.txt.gz",sep="\t", fill = TRUE)## Warning in scan(file = file, what = what, sep = sep, quote = quote, dec = dec, :
## EOF within quoted stringcolnames(go) <- go[1,]
go <- go[2:nrow(go),]
go <- go[which(go$`GO term accession`!=""),]
go_genes <- unique(go$`Gene name`)
length(unique(go$`GO term accession`))## [1] 12264length(go_genes)## [1] 8412v1 <- list("GO"=go_genes, "protein coding"=prot, "non-protein coding"=nprot)
plot(euler(v1),quantities = TRUE)
length(unique(go$`GO term accession`))## [1] 12264head(go)## Gene stable ID Gene name
## 7 ENSG00000198888 MT-ND1
## 8 ENSG00000198888 MT-ND1
## 9 ENSG00000198888 MT-ND1
## 10 ENSG00000198888 MT-ND1
## 11 ENSG00000198888 MT-ND1
## 12 ENSG00000198888 MT-ND1
## GO term definition
## 7 A lipid bilayer along with all the proteins and protein complexes embedded in it an attached to it.
## 8 A metabolic process that results in the removal or addition of one or more electrons to or from a substance, with or without the concomitant removal or addition of a proton or protons.
## 9 The component of a membrane consisting of the gene products and protein complexes having at least some part of their peptide sequence embedded in the hydrophobic region of the membrane.
## 10 The inner, i.e. lumen-facing, lipid bilayer of the mitochondrial envelope. It is highly folded to form cristae.
## 11 A semiautonomous, self replicating organelle that occurs in varying numbers, shapes, and sizes in the cytoplasm of virtually all eukaryotic cells. It is notably the site of tissue respiration.
## 12 Interacting selectively and non-covalently with any protein or protein complex (a complex of two or more proteins that may include other nonprotein molecules).
## GO term name GO term accession GO domain
## 7 membrane GO:0016020 cellular_component
## 8 oxidation-reduction process GO:0055114 biological_process
## 9 integral component of membrane GO:0016021 cellular_component
## 10 mitochondrial inner membrane GO:0005743 cellular_component
## 11 mitochondrion GO:0005739 cellular_component
## 12 protein binding GO:0005515 molecular_function
## GO term evidence code
## 7 IEA
## 8 IEA
## 9 IEA
## 10 IEA
## 11 IEA
## 12 IPIHuman Phenotype Ontology
Now same for Human phenotype ontology. Downloaded from https://hpo.jax.org/app/download/annotation. There are 2958 sets and 4366 genes as of 27/Sep/2020.
hpo <- read.table("phenotype_to_genes.txt",sep="\t")
hpo_genes <- unique(hpo$V4)
head(hpo)## V1 V2 V3 V4 V5 V6 V7
## 1 HP:0000002 Abnormality of body height 7291 TWIST1 orphadata ORPHA:794
## 2 HP:0000002 Abnormality of body height 50485 SMARCAL1 - mim2gene OMIM:242900
## 3 HP:0000002 Abnormality of body height 1487 CTBP1 - mim2gene OMIM:194190
## 4 HP:0000002 Abnormality of body height 57187 THOC2 - mim2gene OMIM:300957
## 5 HP:0000002 Abnormality of body height 10806 SDCCAG8 orphadata ORPHA:3156
## 6 HP:0000002 Abnormality of body height 84295 PHF6 - mim2gene OMIM:301900length(unique(hpo$V1))## [1] 2958length(hpo_genes)## [1] 4366v1 <- list("HPO"=hpo_genes, "protein coding"=prot, "non-protein coding"=nprot)
plot(euler(v1),quantities = TRUE)
MSigDB
Now for MSigDB version 7.2 accessed 27/Sep/2020. There are 31120 sets and 40044 genes.
msigdb <- gmt_import("msigdb.v7.2.symbols.gmt")
msigdb_genes <- unique(unname(unlist(msigdb)))
length(msigdb)## [1] 31120length(msigdb_genes)## [1] 40044v1 <- list("MSigDB"=msigdb_genes,"protein coding"=prot, "non-protein coding"=nprot)
plot(euler(v1),quantities = TRUE)
Defining the overrepresentation analysis function
Here I'm using the triwise::testEnrichment function which uses a Fisher test under the hood.
Reference van de Laar L, Saelens W, De Prijck S, Martens L, Scott CL, Van Isterdael G, Hoffmann E, Beyaert R, Saeys Y, Lambrecht BN, Guilliams M. Yolk Sac Macrophages, Fetal Liver, and Adult Monocytes Can Colonize an Empty Niche and Develop into Functional Tissue-Resident Macrophages. Immunity. 2016 Apr 19;44(4):755-68. doi: 10.1016/j.immuni.2016.02.017. Epub 2016 Mar 15. PMID: 26992565.
ora <- function(gs,genesets,universe){
res <- triwise::testEnrichment(gs, genesets, universe, minknown = 3, mindiffexp = 2, maxknown = 5000)
if( nrow(res)>0 ) {
res <- res[order(res$pval),]
res <- head(res,5)
hits <- log(res[which(res$qval<0.01),2])
names(hits) <- res[which(res$qval<0.01),4]
return(hits)
} else {
return(0)
}
}Here is a function that converts list of vectors to a network diagram. There are 2x edges than nodes. Only the edges with highest similarity are retained, as per jaccard. The size of the gene set is proportional to the node size (sqrt).
gs2net <- function(gset){
mydf <- bind_rows(lapply(gset, as.data.frame.list))
rownames(mydf) <- names(gset)
j <- apply(mydf,1,function(x) {
apply(mydf,1,function(y) {
length(intersect(x,y) ) / length(union(x,y))
})
})
j[lower.tri(j)] <- NA
j[lower.tri(j,diag=TRUE)] <- 0
jl <- melt(j)
jl <- jl[which(jl$Var1 != jl$Var2),]
jl <- jl[which(jl$value != 1),]
jl <- jl[order(-jl$value),]
jl <- head(jl,length(gset)*2)
jl$edgeSize = with(jl, jl$value/sum(jl$value))
lengths <- unlist(lapply(gset,length))
lengths <- sqrt(lengths/sum(lengths)*100)
jl$vertexsize <- lengths[match(as.character(jl$Var1),names(lengths))]
jlnet = network(jl[,1:2])
plot(jlnet, displaylabels = TRUE, label.col = "steelblue",
edge.lwd = c(jl$edgeSize) * 100,
arrowhead.cex = 0,
label.cex = 1, vertex.border = "white",vertex.cex = jl$vertexsize,
vertex.col = "blue", edge.col = rgb(0, 0, 0, alpha = 0.5))
}Epilepsy
Let's analyse the new epilepsy gene sets. There are only 12 sets. We begin with a histogram of the set size.
epi <- gmt_import("epilepsy_genesymbols.gmt")
epi <- diab[which(!names(epi)=='')]
epi <- epi[which(!duplicated(names(epi)))]
length(epi)## [1] 12setnames <- names(epi)
names(epi) <-paste("GS",1:length(epi))
numgenes <- unlist(lapply(epi,length))
hist(numgenes,breaks = 15, xlab="number of genes per set",main="Epilepsy gene sets")
Next we can look at the similarity between the gene sets with a network diagram. GS8, GS12 and GS10 have some overlap, as do GS7 with GS11.
Next we can see that the gene sets have diverse representation of biotypes. There is a high proportion of protein coding genes but also lincRNA.
gs2net(epi)
# classes
mytable <- lapply(epi,function(x) {
table(biotypes_df[which(x %in% biotypes_df$universe),2])
})
mydf <- bind_rows(lapply(mytable, as.data.frame.list))
rownames(mydf) <- names(epi)
mydf_n <- mydf/rowSums(mydf)
mydf <- mydf[,head(order(-colSums(mydf_n)),8)]
colfunc <- colorRampPalette(c("white", "yellow","orange","red","darkred"))(n=25)
heatmap.2(t(as.matrix(mydf)),col=colfunc,trace="none",margin=c(5,10),cexCol = 0.9, cexRow = 0.9,
scale="col",dendrogram = "none",main="Epilepsy gene sets by biotype",na.color = "white")
The overall biotype has a high representaiton of protein coding genes similar to Reactome. We also observed a higher representation of lincRNAs in these gene sets. Pseudogenes were underrepresented.
mydf_n[is.na(mydf_n)] <- 0
mytable <-colSums(mydf_n)
mytable <- mytable/sum(mytable)*100
mytable<-mytable[order(mytable)]
par(mar=c(3,13,1,1)); barplot(mytable,horiz=TRUE,las=2,cex.names = 0.7,cex.axis = 0.8,
main="Overall biotype class representation (%)" ) ;grid()
Next we wanted to know the Reactome functional categories that are enriched in each of the gene sets. Using a Fisher text, we performed an enrichment analysis with Reactome genes using Ensembl as the background (universe). We are aware that it may the improper background and the best background would be all genes that were detected in the original experiment. As expected, some of these sets are related to neurological functions.
# ORA
res <- lapply(epi, function(x) { ora(gs=x,genesets=reactome, universe=universe) } )
res <- res[which(lapply(res,length)>0)]
myres <- bind_rows(lapply(res, as.data.frame.list))
myres[is.na(myres)] <- 0
rownames(myres) <- names(res)
colfunc <- colorRampPalette(c("white", "lightblue","blue", "darkblue"))(n=25)
heatmap.2(t(as.matrix(myres)),col=colfunc,trace="none",margin=c(10,25),cexRow = 0.6,
dendrogram = "none",main="Epilepsy gene sets by Reactome enrichment")
Diabetes
Now lets analyse the new diabetes gene sets. There are 129 gene sets.
diab <- gmt_import("diabetes_genesymbols.gmt")## Warning in gmt_import("diabetes_genesymbols.gmt"): Duplicated gene sets names
## detecteddiab <- diab[which(!names(diab)=='')]
diab <- diab[which(!duplicated(names(diab)))]
length(diab)## [1] 124setnames <- names(diab)
names(diab) <-paste("GS",1:length(diab))
numgenes <- unlist(lapply(diab,length))
hist(numgenes,breaks = 15, xlab="number of genes per set",main="Diabetes gene sets")
The network diagram indicates clustering into a few groups. A large group with GS39 at the centre (right side). A smaller cluster around GS121 (left side). A smaller cluster around GS120 (top).
The biotype heatmap indicates representation of protein coding, lincRNAs and processed pseudgenes. In some sets, non-coding genes outnumber protein coding genes (left). In most sets, protein coding genes are more abundant (right).
gs2net(diab)
# classes
mytable <- lapply(diab,function(x) {
table(biotypes_df[which(x %in% biotypes_df$universe),2])
})
mydf <- bind_rows(lapply(mytable, as.data.frame.list))
rownames(mydf) <- names(diab)
mydf_n <- mydf/rowSums(mydf)
mydf <- mydf[,head(order(-colSums(mydf_n)),8)]
colfunc <- colorRampPalette(c("white", "yellow","orange","red","darkred"))(n=25)
heatmap.2(t(as.matrix(mydf)),col=colfunc,trace="none",margin=c(5,10),cexCol = 0.3,
cexRow = 0.6 ,scale="col",dendrogram = "none",main="Diabetes gene sets by biotype",na.color = "white")
In these gene sets, psudogenes are underrepresented while protein coding, lincRNA and antisense RNA genes are overrepresented.
mydf_n[is.na(mydf_n)] <- 0
mytable <-colSums(mydf_n)
mytable <- mytable/sum(mytable)*100
mytable<-mytable[order(mytable)]
par(mar=c(3,13,1,1)); barplot(mytable,horiz=TRUE,las=2,cex.names = 0.7,cex.axis = 0.8,
main="Overall biotype class representation (%)" ) ;grid()
Here is the reactome enrichment plot for the diabetes gene sets. The text is too small to make out from this chart so I will need to look at making a PDF version. One thing we can see from the chart is that not many of the enriched Reactomes are common apart from the horizontal strip near the bottom.
# reactome enrichment
res <- lapply(diab, function(x) { ora(gs=x,genesets=reactome, universe=universe) } )
res <- res[which(lapply(res,length)>0)]
myres <- bind_rows(lapply(res, as.data.frame.list))
myres[is.na(myres)] <- 0
rownames(myres) <- names(res)
colfunc <- colorRampPalette(c("white", "lightblue","blue", "darkblue"))(n=25)
heatmap.2(t(as.matrix(myres)),col=colfunc,trace="none",margin=c(5,10),cexRow = 0.2,
cexCol = 0.3 ,dendrogram = "none",main="Diabetes gene sets by Reactome enrichment")
Heart disease
Let's analyse the new heart disease gene sets. There are 50 gene sets.
hd <- gmt_import("heartdisease_genesymbols.gmt")
hd <- diab[which(!names(hd)=='')]
hd <- hd[which(!duplicated(names(hd)))]
length(hd)## [1] 50setnames <- names(hd)
names(hd) <-paste("GS",1:length(hd))
numgenes <- unlist(lapply(hd,length))
hist(numgenes,breaks = 15, xlab="number of genes per set",main="Heart disease gene sets")
The network chart indicatesmost sets are linked around GS50. There are several very small groups on the right of the chart.
The biotype heatmap shows overall genes sets are composed motly of protein coding genes and a small number of sets have a majority of lincRNA and pseudo genes.
gs2net(hd)
# classes
mytable <- lapply(hd,function(x) {
table(biotypes_df[which(x %in% biotypes_df$universe),2])
})
mydf <- bind_rows(lapply(mytable, as.data.frame.list))
rownames(mydf) <- names(hd)
mydf_n <- mydf/rowSums(mydf)
mydf <- mydf[,head(order(-colSums(mydf_n)),8)]
colfunc <- colorRampPalette(c("white", "yellow","orange","red","darkred"))(n=25)
heatmap.2(t(as.matrix(mydf)),col=colfunc,trace="none",margin=c(5,10),cexCol = 0.2,
cexRow = 0.8,scale="col",dendrogram = "none",main="Heart disease gene sets by biotype",na.color = "white")
Similar to the other sets, pseudogenes are underrepresented and protein coding, lincRNA, and antisense genes are overrepresented.
mydf_n[is.na(mydf_n)] <- 0
mytable <-colSums(mydf_n)
mytable <- mytable/sum(mytable)*100
mytable<-mytable[order(mytable)]
par(mar=c(3,13,1,1)); barplot(mytable,horiz=TRUE,las=2,cex.names = 0.7,cex.axis = 0.8,
main="Overall biotype class representation (%)" ) ;grid()
The reactome heatmap shows a cluster gene sets that are enriched for ECM, collagen and IGF signaling (lower). There is a strip showing metabolism, immune signaling and signal transduction are common. Apart from that, the enriched gene sets are quite distinct.
# reactome enrichment
res <- lapply(hd, function(x) { ora(gs=x,genesets=reactome, universe=universe) } )
res <- res[which(lapply(res,length)>0)]
myres <- bind_rows(lapply(res, as.data.frame.list))
myres[is.na(myres)] <- 0
rownames(myres) <- names(res)
colfunc <- colorRampPalette(c("white", "lightblue","blue", "darkblue"))(n=25)
heatmap.2(t(as.matrix(myres)),col=colfunc,trace="none",margin=c(5,15),cexRow = 0.4, cexCol = 0.5,
dendrogram = "none",main="Heart disease gene sets by Reactome enrichment")
SARS, MERS and SARS-CoV-2
Let's analyse the new virus infection gene sets. There are 66 sets.
sars <- gmt_import("sarsmers_genesymbols.gmt")## Warning in gmt_import("sarsmers_genesymbols.gmt"): Duplicated gene sets names
## detectedsars <- diab[which(!names(sars)=='')]
sars <- sars[which(!duplicated(names(sars)))]
length(sars)## [1] 65setnames <- names(sars)
names(sars) <-paste("GS",1:length(sars))
numgenes <- unlist(lapply(sars,length))
hist(numgenes,breaks = 15, xlab="number of genes per set",main="SARS, MERS and SARS-CoV-2 gene sets")
The network diagram shows a large cluster of sets around GS21 (left side). A smaller cluster exists around GS24 (right side) and GSGS29 (bottom). The biotype chart shows most gene sets are dominated with protein coding genes but about 1/3 have dominance of lincRNA and pseudogenes.
gs2net(sars)
# classes
mytable <- lapply(sars,function(x) {
table(biotypes_df[which(x %in% biotypes_df$universe),2])
})
mydf <- bind_rows(lapply(mytable, as.data.frame.list))
rownames(mydf) <- names(sars)
mydf_n <- mydf/rowSums(mydf)
mydf <- mydf[,head(order(-colSums(mydf_n)),8)]
colfunc <- colorRampPalette(c("white", "yellow","orange","red","darkred"))(n=25)
heatmap.2(t(as.matrix(mydf)),col=colfunc,trace="none",margin=c(5,10),cexCol = 0.4, cexRow = 0.8,
scale="col",dendrogram = "none",main="SARS, MERS and SARS-CoV2 gene sets by biotype",na.color = "white")
The overall biotype representation is similar to the other sets with overrepresentation of lincRNA and antisense RNA genes.
mydf_n[is.na(mydf_n)] <- 0
mytable <-colSums(mydf_n)
mytable <- mytable/sum(mytable)*100
mytable<-mytable[order(mytable)]
par(mar=c(3,13,1,1)); barplot(mytable,horiz=TRUE,las=2,cex.names = 0.7,cex.axis = 0.8,
main="Overall biotype class representation (%)" ) ;grid()
There is a cluster of immune related enrichment in several gene sets as expected (lower left). Moreover there are some clustering of integrins (lower left), and cell cycle (middle). Post translational modification and signal transduction were also common.
# Reactome enrichment
res <- lapply(sars, function(x) { ora(gs=x,genesets=reactome, universe=universe) } )
res <- res[which(lapply(res,length)>0)]
myres <- bind_rows(lapply(res, as.data.frame.list))
myres[is.na(myres)] <- 0
rownames(myres) <- names(res)
colfunc <- colorRampPalette(c("white", "lightblue","blue", "darkblue"))(n=25)
heatmap.2(t(as.matrix(myres)),col=colfunc,trace="none",margin=c(5,10),cexRow = 0.3, cexCol = 0.4,
dendrogram = "none",main="SARS, MERS and SARS-CoV-2 gene sets genes by Reactome enrichment")
Overall analysis
There are 257 gene sets in total. The median set size was smallest for epilepsy and largest for diabetes. In terms of number of genes in each library, epilepsy has 2175, diabetes has 11179, heart disease has 8727 and virus has 6827. There were 731 genes that were common to all three libraries. Altogether these new sets include 15522 genes, of these 12824 are protein coding and 2698 are non-protein coding.
names(epi) <- paste("Epilepsy",names(epi) )
names(diab) <- paste("Diabetes",names(diab) )
names(hd) <- paste("Heart disease",names(hd) )
names(sars) <- paste("Virus related",names(sars) )
all <- c(epi,diab,hd,sars)
length(all)## [1] 251barplot( c(length(epi) , length(diab) , length(hd) , length(sars) ) ,
names.arg = c("Epilepsy","Diabetes","Heart Disease","Virus") ,
ylab= "number of gene sets")
epil <- sapply(epi,length)
diabl <- sapply(diab,length)
hdl <- sapply(hd,length)
sarsl <- sapply(sars,length)
boxplot(epil,diabl,hdl,sarsl,names=c("Epilepsy","Diabetes","Heart Disease","Virus"))
epig <- unique(unname(unlist(epi)))
diabg <- unique(unname(unlist(diab)))
hdg <- unique(unname(unlist(hd)))
sarsg <- unique(unname(unlist(sars)))
length(epig)## [1] 1373length(diabg)## [1] 11164length(hdg)## [1] 7060length(sarsg)## [1] 7472v1 <- list("Epilepsy"=epig, "Diabetes"=diabg, "Heart diabetes"=hdg, "Virus related" = sarsg)
plot(euler(v1),quantities = TRUE)
allg <- unique(unname(unlist(all)))
length(allg)## [1] 11164v1 <- list("DEE2 sets"=allg,"protein coding"=prot, "non-protein coding"=nprot)
plot(euler(v1),quantities = TRUE)
Session information
sessionInfo()## R version 4.0.2 (2020-06-22)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.5 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1
##
## locale:
## [1] LC_CTYPE=en_AU.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_AU.UTF-8 LC_COLLATE=en_AU.UTF-8
## [5] LC_MONETARY=en_AU.UTF-8 LC_MESSAGES=en_AU.UTF-8
## [7] LC_PAPER=en_AU.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_AU.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] parallel stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] eulerr_6.1.0 network_1.16.1 reshape2_1.4.4
## [4] gplots_3.1.0 dplyr_1.0.2 triwise_0.99.5
## [7] Biobase_2.48.0 BiocGenerics_0.34.0 mitch_1.1.12
## [10] getDEE2_0.99.30
##
## loaded via a namespace (and not attached):
## [1] Rcpp_1.0.5 lattice_0.20-41
## [3] gtools_3.8.2 digest_0.6.25
## [5] mime_0.9 R6_2.4.1
## [7] GenomeInfoDb_1.24.2 plyr_1.8.6
## [9] stats4_4.0.2 evaluate_0.14
## [11] ggplot2_3.3.2 pillar_1.4.6
## [13] zlibbioc_1.34.0 rlang_0.4.8
## [15] S4Vectors_0.26.1 Matrix_1.2-18
## [17] rmarkdown_2.4 stringr_1.4.0
## [19] htmlwidgets_1.5.2 polyclip_1.10-0
## [21] RCurl_1.98-1.2 munsell_0.5.0
## [23] shiny_1.5.0 DelayedArray_0.14.1
## [25] polylabelr_0.2.0 compiler_4.0.2
## [27] httpuv_1.5.4 xfun_0.18
## [29] pkgconfig_2.0.3 htmltools_0.5.0
## [31] tidyselect_1.1.0 SummarizedExperiment_1.18.2
## [33] tibble_3.0.3 gridExtra_2.3
## [35] GenomeInfoDbData_1.2.3 IRanges_2.22.2
## [37] matrixStats_0.57.0 reshape_0.8.8
## [39] crayon_1.3.4 later_1.1.0.1
## [41] MASS_7.3-53 bitops_1.0-6
## [43] grid_4.0.2 xtable_1.8-4
## [45] GGally_2.0.0 gtable_0.3.0
## [47] lifecycle_0.2.0 magrittr_1.5
## [49] scales_1.1.1 KernSmooth_2.23-17
## [51] stringi_1.5.3 XVector_0.28.0
## [53] promises_1.1.1 ellipsis_0.3.1
## [55] generics_0.0.2 vctrs_0.3.4
## [57] htm2txt_2.1.1 pbmcapply_1.5.0
## [59] RColorBrewer_1.1-2 tools_4.0.2
## [61] glue_1.4.2 beeswarm_0.2.3
## [63] purrr_0.3.4 yaml_2.2.1
## [65] fastmap_1.0.1 colorspace_1.4-1
## [67] GenomicRanges_1.40.0 caTools_1.18.0
## [69] knitr_1.30 echarts4r_0.3.2