BiocSet 1.2.2
BiocSet is a package that represents element sets in a tibble format with the
BiocSet class. Element sets are read in and converted into a tibble format.
From here, typical dplyr operations can be performed on the element set.
BiocSet also provides functionality for mapping different ID types and
providing reference urls for elements and sets.
Install the most recent version from Bioconductor:
if(!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("BiocSet")The development version is also available for install from GitHub:
BiocManager::install("Kayla-Morrell/BiocSet")Then load BiocSet:
library(BiocSet)BiocSet can create a BiocSet object using two different input methods. The
first is to input named character vectors of element sets. BiocSet returns
three tibbles, es_element which contains the elements, es_set which contains
the sets and es_elementset which contains elements and sets together.
tbl <- BiocSet(set1 = letters, set2 = LETTERS)
tbl
#> class: BiocSet
#>
#> es_element():
#> # A tibble: 52 x 1
#> element
#> <chr>
#> 1 a
#> 2 b
#> 3 c
#> # ... with 49 more rows
#>
#> es_set():
#> # A tibble: 2 x 1
#> set
#> <chr>
#> 1 set1
#> 2 set2
#>
#> es_elementset() <active>:
#> # A tibble: 52 x 2
#> element set
#> <chr> <chr>
#> 1 a set1
#> 2 b set1
#> 3 c set1
#> # ... with 49 more rowsThe second method of creating a BiocSet object would be to read in a .gmt
file. Using import(), a path to a downloaded .gmt file is read in and a
BiocSet object is returned. The example below uses a hallmark element set
downloaded from GSEA, which is also included with this package. This
BiocSet includes a source column within the es_elementset tibble for
reference as to where the element set came from.
gmtFile <- system.file(package = "BiocSet",
"extdata",
"hallmark.gene.symbol.gmt")
tbl2 <- import(gmtFile)
tbl2
#> class: BiocSet
#>
#> es_element():
#> # A tibble: 4,386 x 1
#> element
#> <chr>
#> 1 JUNB
#> 2 CXCL2
#> 3 ATF3
#> # ... with 4,383 more rows
#>
#> es_set():
#> # A tibble: 50 x 2
#> set source
#> <chr> <chr>
#> 1 HALLMARK_TNFA_SIGNALIN~ http://www.broadinstitute.org/gsea/msigdb/cards/HALLM~
#> 2 HALLMARK_HYPOXIA http://www.broadinstitute.org/gsea/msigdb/cards/HALLM~
#> 3 HALLMARK_CHOLESTEROL_H~ http://www.broadinstitute.org/gsea/msigdb/cards/HALLM~
#> # ... with 47 more rows
#>
#> es_elementset() <active>:
#> # A tibble: 7,324 x 2
#> element set
#> <chr> <chr>
#> 1 JUNB HALLMARK_TNFA_SIGNALING_VIA_NFKB
#> 2 CXCL2 HALLMARK_TNFA_SIGNALING_VIA_NFKB
#> 3 ATF3 HALLMARK_TNFA_SIGNALING_VIA_NFKB
#> # ... with 7,321 more rowsexport() allows for a BiocSet object to be exported into a temporary file
with the extention .gmt.
fl <- tempfile(fileext = ".gmt")
gmt <- export(tbl2, fl)
gmt
#> GMTFile object
#> resource: C:\Users\biocbuild\bbs-3.11-bioc\tmpdir\RtmpimjyXI\file2c24ea44469.gmtA feature available to BiocSet is the ability to activate different
tibbles to perform certain functions on. When a BiocSet is created, the tibble
es_elementset is automatically activated and all functions will be performed
on this tibble. BiocSet adopts the use of many common dplyr functions such
as filter(), select(), mutate(), summarise(), and arrange(). With
each of the functions the user is able to pick a different tibble to activate
and work on by using ‘verb_tibble’. After the function is executed than the
‘active’ tibble is returned back to the tibble that was active before the
function call. Some examples are shown below of how these functions work.
tbl <- BiocSet(set1 = letters, set2 = LETTERS)
tbl
#> class: BiocSet
#>
#> es_element():
#> # A tibble: 52 x 1
#> element
#> <chr>
#> 1 a
#> 2 b
#> 3 c
#> # ... with 49 more rows
#>
#> es_set():
#> # A tibble: 2 x 1
#> set
#> <chr>
#> 1 set1
#> 2 set2
#>
#> es_elementset() <active>:
#> # A tibble: 52 x 2
#> element set
#> <chr> <chr>
#> 1 a set1
#> 2 b set1
#> 3 c set1
#> # ... with 49 more rows
tbl %>% filter_element(element == "a" | element == "A")
#> class: BiocSet
#>
#> es_element():
#> # A tibble: 2 x 1
#> element
#> <chr>
#> 1 a
#> 2 A
#>
#> es_set():
#> # A tibble: 2 x 1
#> set
#> <chr>
#> 1 set1
#> 2 set2
#>
#> es_elementset() <active>:
#> # A tibble: 2 x 2
#> element set
#> <chr> <chr>
#> 1 a set1
#> 2 A set2
tbl %>% mutate_set(pval = rnorm(1:2))
#> class: BiocSet
#>
#> es_element():
#> # A tibble: 52 x 1
#> element
#> <chr>
#> 1 a
#> 2 b
#> 3 c
#> # ... with 49 more rows
#>
#> es_set():
#> # A tibble: 2 x 2
#> set pval
#> <chr> <dbl>
#> 1 set1 -0.577
#> 2 set2 1.39
#>
#> es_elementset() <active>:
#> # A tibble: 52 x 2
#> element set
#> <chr> <chr>
#> 1 a set1
#> 2 b set1
#> 3 c set1
#> # ... with 49 more rows
tbl %>% arrange_elementset(desc(element))
#> class: BiocSet
#>
#> es_element():
#> # A tibble: 52 x 1
#> element
#> <chr>
#> 1 a
#> 2 b
#> 3 c
#> # ... with 49 more rows
#>
#> es_set():
#> # A tibble: 2 x 1
#> set
#> <chr>
#> 1 set1
#> 2 set2
#>
#> es_elementset() <active>:
#> # A tibble: 52 x 2
#> element set
#> <chr> <chr>
#> 1 z set1
#> 2 y set1
#> 3 x set1
#> # ... with 49 more rowsBiocSet also allows for common set operations to be performed on the BiocSet
object. union() and intersection() are the two set operations available in
BiocSet. We demonstate how a user can find the union between two BiocSet
objects or within a single BiocSet object. Intersection is used in the same
way.
# union of two BiocSet objects
es1 <- BiocSet(set1 = letters[c(1:3)], set2 = LETTERS[c(1:3)])
es2 <- BiocSet(set1 = letters[c(2:4)], set2 = LETTERS[c(2:4)])
union(es1, es2)
#> class: BiocSet
#>
#> es_element():
#> # A tibble: 8 x 1
#> element
#> <chr>
#> 1 a
#> 2 b
#> 3 c
#> # ... with 5 more rows
#>
#> es_set():
#> # A tibble: 2 x 1
#> set
#> <chr>
#> 1 set1
#> 2 set2
#>
#> es_elementset() <active>:
#> # A tibble: 8 x 2
#> element set
#> <chr> <chr>
#> 1 a set1
#> 2 b set1
#> 3 c set1
#> # ... with 5 more rows
# union within a single BiocSet object
es3 <- BiocSet(set1 = letters[c(1:10)], set2 = letters[c(4:20)])
union_single(es3)
#> class: BiocSet
#>
#> es_element():
#> # A tibble: 20 x 1
#> element
#> <chr>
#> 1 a
#> 2 b
#> 3 c
#> # ... with 17 more rows
#>
#> es_set():
#> # A tibble: 1 x 1
#> set
#> <chr>
#> 1 union
#>
#> es_elementset() <active>:
#> # A tibble: 20 x 2
#> element set
#> <chr> <chr>
#> 1 a union
#> 2 b union
#> 3 c union
#> # ... with 17 more rowsNext, we demonstrate the a couple uses of BiocSet with an experiment dataset
airway from the package airway. This data is from an RNA-Seq experiment on
airway smooth muscle (ASM) cell lines.
The first step is to load the library and the necessary data.
library(airway)
data("airway")
se <- airwayThe function go_sets() discovers the keys from the org object and uses
AnnotationDbi::select to create a mapping to GO ids. go_sets() also allows
the user to indicate which evidence type or ontology type they would like when
selecting the GO ids. The default is using all evidence types and all ontology
types. We represent these identifieres as a BiocSet object. Using the
go_sets function we are able to map the Ensembl ids and GO ids from the genome
wide annotation for Human data in the org.Hs.eg.db package. The Ensembl ids
are treated as elements while the GO ids are treated as sets.
library(org.Hs.eg.db)
go <- go_sets(org.Hs.eg.db, "ENSEMBL")
go
#> class: BiocSet
#>
#> es_element():
#> # A tibble: 22,432 x 1
#> element
#> <chr>
#> 1 ENSG00000151729
#> 2 ENSG00000025708
#> 3 ENSG00000068305
#> # ... with 22,429 more rows
#>
#> es_set():
#> # A tibble: 18,105 x 1
#> set
#> <chr>
#> 1 GO:0000002
#> 2 GO:0000003
#> 3 GO:0000010
#> # ... with 18,102 more rows
#>
#> es_elementset() <active>:
#> # A tibble: 317,035 x 2
#> element set
#> <chr> <chr>
#> 1 ENSG00000151729 GO:0000002
#> 2 ENSG00000025708 GO:0000002
#> 3 ENSG00000068305 GO:0000002
#> # ... with 317,032 more rows
# an example of subsetting by evidence type
go_sets(org.Hs.eg.db, "ENSEMBL", evidence = c("IPI", "TAS"))
#> class: BiocSet
#>
#> es_element():
#> # A tibble: 15,976 x 1
#> element
#> <chr>
#> 1 ENSG00000151729
#> 2 ENSG00000168685
#> 3 ENSG00000114030
#> # ... with 15,973 more rows
#>
#> es_set():
#> # A tibble: 4,305 x 2
#> set evidence
#> <chr> <named list>
#> 1 GO:0000002 <chr [1]>
#> 2 GO:0000018 <chr [1]>
#> 3 GO:0000019 <chr [1]>
#> # ... with 4,302 more rows
#>
#> es_elementset() <active>:
#> # A tibble: 67,915 x 2
#> element set
#> <chr> <chr>
#> 1 ENSG00000151729 GO:0000002
#> 2 ENSG00000168685 GO:0000018
#> 3 ENSG00000114030 GO:0000018
#> # ... with 67,912 more rowsSome users may not be interested in reporting the non-descriptive elements. We
demonstrate subsetting the airway data to include non-zero assays and then
filter out the non-descriptive elements.
se1 = se[rowSums(assay(se)) != 0,]
go %>% filter_element(element %in% rownames(se1))
#> class: BiocSet
#>
#> es_element():
#> # A tibble: 16,472 x 1
#> element
#> <chr>
#> 1 ENSG00000151729
#> 2 ENSG00000025708
#> 3 ENSG00000068305
#> # ... with 16,469 more rows
#>
#> es_set():
#> # A tibble: 18,105 x 1
#> set
#> <chr>
#> 1 GO:0000002
#> 2 GO:0000003
#> 3 GO:0000010
#> # ... with 18,102 more rows
#>
#> es_elementset() <active>:
#> # A tibble: 249,545 x 2
#> element set
#> <chr> <chr>
#> 1 ENSG00000151729 GO:0000002
#> 2 ENSG00000025708 GO:0000002
#> 3 ENSG00000068305 GO:0000002
#> # ... with 249,542 more rowsIt may also be of interest to users to know how many elements are in each set.
Using the count function we are able to calculate the elements per set.
go %>% group_by(set) %>% dplyr::count()
#> Warning: The `add` argument of `group_by()` is deprecated as of dplyr 1.0.0.
#> Please use the `.add` argument instead.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_warnings()` to see where this warning was generated.
#> # A tibble: 18,105 x 2
#> # Groups: set [18,105]
#> set n
#> <chr> <int>
#> 1 GO:0000002 13
#> 2 GO:0000003 4
#> 3 GO:0000010 2
#> 4 GO:0000012 7
#> 5 GO:0000014 9
#> 6 GO:0000015 4
#> 7 GO:0000016 1
#> 8 GO:0000018 6
#> 9 GO:0000019 3
#> 10 GO:0000022 2
#> # ... with 18,095 more rowsIt may also be helpful to remove sets that are empty. Since we have shown how to calculate the number of elements per set, we know that this data set does not contain any empty sets. We decide to demonstrate regardless for those users that may need this functionality.
drop <- es_activate(go, elementset) %>% group_by(set) %>%
dplyr::count() %>% filter(n == 0) %>% pull(set)
go %>% filter_set(!(set %in% drop))
#> class: BiocSet
#>
#> es_element():
#> # A tibble: 22,432 x 1
#> element
#> <chr>
#> 1 ENSG00000151729
#> 2 ENSG00000025708
#> 3 ENSG00000068305
#> # ... with 22,429 more rows
#>
#> es_set():
#> # A tibble: 18,105 x 1
#> set
#> <chr>
#> 1 GO:0000002
#> 2 GO:0000003
#> 3 GO:0000010
#> # ... with 18,102 more rows
#>
#> es_elementset() <active>:
#> # A tibble: 317,035 x 2
#> element set
#> <chr> <chr>
#> 1 ENSG00000151729 GO:0000002
#> 2 ENSG00000025708 GO:0000002
#> 3 ENSG00000068305 GO:0000002
#> # ... with 317,032 more rowsTo simplify mapping we created a couple map functions. map_unique() is used
when there is known 1:1 mapping, This takes four arguements, a BiocSet
object, an AnnotationDbi object, the id to map from, and the id to map to.
map_multiple() needs a fifth argument to indicate how the function should
treat an element when there is multiple mapping. Both functions utilize
mapIds from AnnotationDbi and return a BiocSet object. In the example
below we show how to use map_unique to map go’s ids from Ensembl to gene
symbols.
go %>% map_unique(org.Hs.eg.db, "ENSEMBL", "SYMBOL")
#> 'select()' returned 1:many mapping between keys and columns
#> Joining, by = "element"
#> `mutate_if()` ignored the following grouping variables:
#> Column `element`
#> Joining, by = "element"
#> Joining, by = "element"
#> Joining, by = "set"
#> Joining, by = c("element", "set")
#> class: BiocSet
#>
#> es_element():
#> # A tibble: 19,821 x 1
#> element
#> <chr>
#> 1 AKT3
#> 2 LONP1
#> 3 MEF2A
#> # ... with 19,818 more rows
#>
#> es_set():
#> # A tibble: 18,105 x 1
#> set
#> <chr>
#> 1 GO:0000002
#> 2 GO:0000003
#> 3 GO:0000010
#> # ... with 18,102 more rows
#>
#> es_elementset() <active>:
#> # A tibble: 280,881 x 2
#> element set
#> <chr> <chr>
#> 1 AKT3 GO:0000002
#> 2 LONP1 GO:0000002
#> 3 MEF2A GO:0000002
#> # ... with 280,878 more rowsAnother functionality of BiocSet is the ability to add information to the
tibbles. Using the GO.db library we are able to map definitions to the GO
ids. From there we can add the mapping to the tibble using map_add() and the
mutate function.
library(GO.db)
map <- map_add_set(go, GO.db, "GOID", "DEFINITION")
go %>% mutate_set(definition = map)
#> class: BiocSet
#>
#> es_element():
#> # A tibble: 22,432 x 1
#> element
#> <chr>
#> 1 ENSG00000151729
#> 2 ENSG00000025708
#> 3 ENSG00000068305
#> # ... with 22,429 more rows
#>
#> es_set():
#> # A tibble: 18,105 x 2
#> set definition
#> <chr> <chr>
#> 1 GO:0000002 The maintenance of the structure and integrity of the mitochondria~
#> 2 GO:0000003 The production of new individuals that contain some portion of gen~
#> 3 GO:0000010 Catalysis of the reaction: all-trans-hexaprenyl diphosphate + isop~
#> # ... with 18,102 more rows
#>
#> es_elementset() <active>:
#> # A tibble: 317,035 x 2
#> element set
#> <chr> <chr>
#> 1 ENSG00000151729 GO:0000002
#> 2 ENSG00000025708 GO:0000002
#> 3 ENSG00000068305 GO:0000002
#> # ... with 317,032 more rowsThe library KEGGREST is a client interface to the KEGG REST server.
KEGG contains pathway maps that represent interaction, reaction and relation
networks for various biological processes and diseases. BiocSet has a
function that utilizes KEGGREST to develop a BiocSet object that contains
the elements for every pathway map in KEGG.
Due to limiations of the KEGGREST package, kegg_sets can take some time to
run depending on the amount of pathways for the species of interest. Therefore
we demonstrate using BiocFileCache to make the data available to the user.
library(BiocFileCache)
rname <- "kegg_hsa"
exists <- NROW(bfcquery(query=rname, field="rname")) != 0L
if (!exists)
{
kegg <- kegg_sets("hsa")
fl <- bfcnew(rname = rname, ext = ".gmt")
export(kegg_sets("hsa"), fl)
}
kegg <- import(bfcrpath(rname=rname))Within the kegg_sets() function we remove pathways that do not contain any
elements. We then mutate the element tibble using the map_add function to
contain both Ensembl and Entrez ids.
map <- map_add_element(kegg, org.Hs.eg.db, "ENTREZID", "ENSEMBL")
#> 'select()' returned 1:many mapping between keys and columns
kegg <- kegg %>% mutate_element(ensembl = map)Since we are working with ASM data we thought we would subset the airway data
to contain only the elements in the asthma pathway. This filter is performed
on the KEGG id, which for asthma is “hsa05310”.
asthma <- kegg %>% filter_set(set == "hsa05310")
se <- se[rownames(se) %in% es_element(asthma)$ensembl,]
se
#> class: RangedSummarizedExperiment
#> dim: 7704 8
#> metadata(1): ''
#> assays(1): counts
#> rownames(7704): ENSG00000000419 ENSG00000000938 ... ENSG00000273079
#> ENSG00000273085
#> rowData names(0):
#> colnames(8): SRR1039508 SRR1039509 ... SRR1039520 SRR1039521
#> colData names(9): SampleName cell ... Sample BioSampleThe filtering can also be done for multiple pathways.
pathways <- c("hsa05310", "hsa04110", "hsa05224", "hsa04970")
multipaths <- kegg %>% filter_set(set %in% pathways)
multipaths
#> class: BiocSet
#>
#> es_element():
#> # A tibble: 7,931 x 2
#> element ensembl
#> <chr> <chr>
#> 1 3101 ENSG00000160883
#> 2 3098 ENSG00000156515
#> 3 3099 ENSG00000159399
#> # ... with 7,928 more rows
#>
#> es_set():
#> # A tibble: 4 x 2
#> set source
#> <chr> <chr>
#> 1 hsa04110 <NA>
#> 2 hsa04970 <NA>
#> 3 hsa05224 <NA>
#> # ... with 1 more row
#>
#> es_elementset() <active>:
#> # A tibble: 392 x 2
#> element set
#> <chr> <chr>
#> 1 595 hsa04110
#> 2 894 hsa04110
#> 3 896 hsa04110
#> # ... with 389 more rowsThis example will start out the same way that Case Study I started, by loading
in the airway dataset, but we will also do some reformating of the data. The
end goal is to be able to perform a Gene Set Enrichment Analysis on the data and
return a BiocSet object of the gene sets.
data("airway")
airway$dex <- relevel(airway$dex, "untrt")Similar to other analyses we perform a differential expression analysis on the
airway data using the library DESeq2 and then store the results into a
tibble.
library(DESeq2)
library(tibble)
des <- DESeqDataSet(airway, design = ~ cell + dex)
des <- DESeq(des)
#> estimating size factors
#> estimating dispersions
#> gene-wise dispersion estimates
#> mean-dispersion relationship
#> final dispersion estimates
#> fitting model and testing
res <- results(des)
tbl <- res %>%
as.data.frame() %>%
as_tibble(rownames = "ENSEMBL") Since we want to use limma::goana() to perform the GSEA we will need to have
ENTREZ identifiers in the data, as well as filter out some NA information.
Later on this will be considered our ‘element’ tibble.
tbl <- tbl %>%
mutate(
ENTREZID = mapIds(
org.Hs.eg.db, ENSEMBL, "ENTREZID", "ENSEMBL"
) %>% unname()
)
#> 'select()' returned 1:many mapping between keys and columns
tbl <- tbl %>% filter(!is.na(padj), !is.na(ENTREZID))
tbl
#> # A tibble: 14,585 x 8
#> ENSEMBL baseMean log2FoldChange lfcSE stat pvalue padj ENTREZID
#> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
#> 1 ENSG0000000~ 709. -0.381 0.101 -3.79 1.52e-4 1.28e-3 7105
#> 2 ENSG0000000~ 520. 0.207 0.112 1.84 6.53e-2 1.97e-1 8813
#> 3 ENSG0000000~ 237. 0.0379 0.143 0.264 7.92e-1 9.11e-1 57147
#> 4 ENSG0000000~ 57.9 -0.0882 0.287 -0.307 7.59e-1 8.95e-1 55732
#> 5 ENSG0000000~ 5817. 0.426 0.0883 4.83 1.38e-6 1.82e-5 3075
#> 6 ENSG0000000~ 1282. -0.241 0.0887 -2.72 6.58e-3 3.28e-2 2519
#> 7 ENSG0000000~ 610. -0.0476 0.167 -0.286 7.75e-1 9.03e-1 2729
#> 8 ENSG0000000~ 369. -0.500 0.121 -4.14 3.48e-5 3.42e-4 4800
#> 9 ENSG0000000~ 183. -0.124 0.180 -0.689 4.91e-1 7.24e-1 90529
#> 10 ENSG0000000~ 2814. -0.0411 0.103 -0.400 6.89e-1 8.57e-1 57185
#> # ... with 14,575 more rowsNow that the data is ready for GSEA we can go ahead and use goana() and make
the results into a tibble. This tibble will be considered our ‘set’ tibble.
library(limma)
#>
#> Attaching package: 'limma'
#> The following object is masked from 'package:DESeq2':
#>
#> plotMA
#> The following object is masked from 'package:BiocGenerics':
#>
#> plotMA
go_ids <- goana(tbl$ENTREZID[tbl$padj < 0.05], tbl$ENTREZID, "Hs") %>%
as.data.frame() %>%
as_tibble(rownames = "GOALL")
go_ids
#> # A tibble: 21,148 x 6
#> GOALL Term Ont N DE P.DE
#> <chr> <chr> <chr> <dbl> <dbl> <dbl>
#> 1 GO:00017~ cell activation BP 944 334 8.47e-12
#> 2 GO:00022~ immune effector process BP 819 256 1.66e- 4
#> 3 GO:00022~ cell activation involved in immune resp~ BP 500 148 2.75e- 2
#> 4 GO:00022~ myeloid leukocyte activation BP 465 149 1.23e- 3
#> 5 GO:00022~ myeloid cell activation involved in imm~ BP 399 118 4.58e- 2
#> 6 GO:00022~ neutrophil activation involved in immun~ BP 361 108 4.05e- 2
#> 7 GO:00023~ leukocyte activation involved in immune~ BP 496 148 2.11e- 2
#> 8 GO:00023~ immune system process BP 1960 650 2.32e-15
#> 9 GO:00024~ leukocyte mediated immunity BP 544 167 4.83e- 3
#> 10 GO:00024~ myeloid leukocyte mediated immunity BP 409 124 2.03e- 2
#> # ... with 21,138 more rowsThe last thing we need to do is create a tibble that we will consider our ‘elementset’ tibble. This tibble will be a mapping of all the elements and sets.
foo <- AnnotationDbi::select(
org.Hs.eg.db,
tbl$ENTREZID,
"GOALL",
"ENTREZID") %>% as_tibble()
#> 'select()' returned many:many mapping between keys and columns
foo <- foo %>% dplyr::select(-EVIDENCEALL) %>% distinct()
foo <- foo %>% filter(ONTOLOGYALL == "BP") %>% dplyr::select(-ONTOLOGYALL)
foo
#> # A tibble: 1,149,541 x 2
#> ENTREZID GOALL
#> <chr> <chr>
#> 1 7105 GO:0001816
#> 2 7105 GO:0001817
#> 3 7105 GO:0001819
#> 4 7105 GO:0002221
#> 5 7105 GO:0002252
#> 6 7105 GO:0002376
#> 7 7105 GO:0002682
#> 8 7105 GO:0002683
#> 9 7105 GO:0002697
#> 10 7105 GO:0002698
#> # ... with 1,149,531 more rowsThe function BiocSet_from_elementset() allows for users to create a BiocSet
object from tibbles. This function is helpful when there is metadata contained
in the tibble that should be in the BiocSet object. For this function to work
properly, the columns that are being joined on need to be named correctly. For
instance, in order to use this function on the tibbles we created we need to
change the column in the ‘element’ tibble to ‘element’, the column in the ‘set’
tibble to ‘set’ and the same will be for the ‘elementset’ tibble. We demonstrate
this below and then create the BiocSet object with the simple function.
foo <- foo %>% dplyr::rename(element = ENTREZID, set = GOALL)
tbl <- tbl %>% dplyr::rename(element = ENTREZID)
go_ids <- go_ids %>% dplyr::rename(set = GOALL)
es <- BiocSet_from_elementset(foo, tbl, go_ids)
#> Joining, by = "element"
#> Joining, by = "set"
#> Joining, by = c("element", "set")
#> more elements in 'element' than in 'elementset'
#> more elements in 'set' than in 'elementset'
es
#> class: BiocSet
#>
#> es_element():
#> # A tibble: 12,270 x 8
#> element ENSEMBL baseMean log2FoldChange lfcSE stat pvalue padj
#> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 3980 ENSG00000005~ 206. -0.393 0.172 -2.29 2.23e- 2 8.66e-2
#> 2 1890 ENSG00000025~ 145. 1.23 0.199 6.18 6.58e-10 1.49e-8
#> 3 4205 ENSG00000068~ 1566. 0.529 0.0919 5.76 8.45e- 9 1.64e-7
#> # ... with 12,267 more rows
#>
#> es_set():
#> # A tibble: 15,000 x 6
#> set Term Ont N DE P.DE
#> <chr> <chr> <chr> <dbl> <dbl> <dbl>
#> 1 GO:0000002 mitochondrial genome maintenance BP 21 6 0.466
#> 2 GO:0000003 reproduction BP 963 309 0.00000334
#> 3 GO:0000012 single strand break repair BP 8 1 0.908
#> # ... with 14,997 more rows
#>
#> es_elementset() <active>:
#> # A tibble: 1,149,541 x 2
#> element set
#> <chr> <chr>
#> 1 3980 GO:0000002
#> 2 1890 GO:0000002
#> 3 4205 GO:0000002
#> # ... with 1,149,538 more rowsFor those users that may need to put this metadata filled BiocSet object back into an object similar to GRanges or SummarizedExperiment, we have created functions that allow for the BiocSet object to be created into a tibble or data.frame.
tibble_from_element(es)
#> Joining, by = "set"
#> Joining, by = "element"
#> # A tibble: 12,262 x 14
#> element set Term Ont N DE P.DE ENSEMBL baseMean log2FoldChange
#> <chr> <lis> <lis> <lis> <lis> <lis> <lis> <list> <list> <list>
#> 1 1 <chr~ <chr~ <chr~ <dbl~ <dbl~ <dbl~ <chr [~ <dbl [3~ <dbl [32]>
#> 2 100 <chr~ <chr~ <chr~ <dbl~ <dbl~ <dbl~ <chr [~ <dbl [4~ <dbl [479]>
#> 3 1000 <chr~ <chr~ <chr~ <dbl~ <dbl~ <dbl~ <chr [~ <dbl [2~ <dbl [255]>
#> 4 10000 <chr~ <chr~ <chr~ <dbl~ <dbl~ <dbl~ <chr [~ <dbl [1~ <dbl [159]>
#> 5 10001 <chr~ <chr~ <chr~ <dbl~ <dbl~ <dbl~ <chr [~ <dbl [7~ <dbl [79]>
#> 6 10003 <chr~ <chr~ <chr~ <dbl~ <dbl~ <dbl~ <chr [~ <dbl [3~ <dbl [31]>
#> 7 10004 <chr~ <chr~ <chr~ <dbl~ <dbl~ <dbl~ <chr [~ <dbl [1~ <dbl [19]>
#> 8 100048~ <chr~ <chr~ <chr~ <dbl~ <dbl~ <dbl~ <chr [~ <dbl [9~ <dbl [95]>
#> 9 10005 <chr~ <chr~ <chr~ <dbl~ <dbl~ <dbl~ <chr [~ <dbl [1~ <dbl [152]>
#> 10 10006 <chr~ <chr~ <chr~ <dbl~ <dbl~ <dbl~ <chr [~ <dbl [1~ <dbl [166]>
#> # ... with 12,252 more rows, and 4 more variables: lfcSE <list>, stat <list>,
#> # pvalue <list>, padj <list>
head(data.frame_from_elementset(es))
#> Joining, by = "set"
#> Joining, by = "element"
#> element set Term Ont N DE P.DE
#> 1 3980 GO:0000002 mitochondrial genome maintenance BP 21 6 0.4661969
#> 2 1890 GO:0000002 mitochondrial genome maintenance BP 21 6 0.4661969
#> 3 4205 GO:0000002 mitochondrial genome maintenance BP 21 6 0.4661969
#> 4 10891 GO:0000002 mitochondrial genome maintenance BP 21 6 0.4661969
#> 5 55186 GO:0000002 mitochondrial genome maintenance BP 21 6 0.4661969
#> 6 4358 GO:0000002 mitochondrial genome maintenance BP 21 6 0.4661969
#> ENSEMBL baseMean log2FoldChange lfcSE stat pvalue
#> 1 ENSG00000005156 205.6627 -0.39272067 0.17180717 -2.285822 2.226466e-02
#> 2 ENSG00000025708 144.6965 1.23104228 0.19933257 6.175821 6.582046e-10
#> 3 ENSG00000068305 1565.5675 0.52905471 0.09186326 5.759155 8.453618e-09
#> 4 ENSG00000109819 152.4021 -0.23932270 0.17393025 -1.375969 1.688311e-01
#> 5 ENSG00000114120 2065.4720 0.14500117 0.09361789 1.548862 1.214150e-01
#> 6 ENSG00000115204 1091.6398 -0.09563418 0.09191349 -1.040480 2.981168e-01
#> padj
#> 1 8.661787e-02
#> 2 1.494473e-08
#> 3 1.636969e-07
#> 4 3.845953e-01
#> 5 3.059216e-01
#> 6 5.485252e-01A final feature to BiocSet is the ability to add reference information about
all of the elements/sets. A user could utilize the function url_ref() to add
information to the BiocSet object. If a user has a question about a specific
id then they can follow the reference url to get more informtion. Below is an
example of using url_ref() to add reference urls to the go data set we
worked with above.
url_ref(go)
#> class: BiocSet
#>
#> es_element():
#> # A tibble: 22,432 x 2
#> element url
#> <chr> <chr>
#> 1 ENSG000001517~ https://www.ensembl.org/Homo_sapiens/Gene/Summary?db=core;g=EN~
#> 2 ENSG000000257~ https://www.ensembl.org/Homo_sapiens/Gene/Summary?db=core;g=EN~
#> 3 ENSG000000683~ https://www.ensembl.org/Homo_sapiens/Gene/Summary?db=core;g=EN~
#> # ... with 22,429 more rows
#>
#> es_set():
#> # A tibble: 18,105 x 2
#> set url
#> <chr> <chr>
#> 1 GO:0000002 http://amigo.geneontology.org/amigo/medial_search?q=GO:0000002
#> 2 GO:0000003 http://amigo.geneontology.org/amigo/medial_search?q=GO:0000003
#> 3 GO:0000010 http://amigo.geneontology.org/amigo/medial_search?q=GO:0000010
#> # ... with 18,102 more rows
#>
#> es_elementset() <active>:
#> # A tibble: 317,035 x 2
#> element set
#> <chr> <chr>
#> 1 ENSG00000151729 GO:0000002
#> 2 ENSG00000025708 GO:0000002
#> 3 ENSG00000068305 GO:0000002
#> # ... with 317,032 more rowssessionInfo()
#> R version 4.0.2 (2020-06-22)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows Server 2012 R2 x64 (build 9600)
#>
#> Matrix products: default
#>
#> locale:
#> [1] LC_COLLATE=C
#> [2] LC_CTYPE=English_United States.1252
#> [3] LC_MONETARY=English_United States.1252
#> [4] LC_NUMERIC=C
#> [5] LC_TIME=English_United States.1252
#>
#> attached base packages:
#> [1] parallel stats4 stats graphics grDevices utils datasets
#> [8] methods base
#>
#> other attached packages:
#> [1] limma_3.44.3 tibble_3.0.3
#> [3] DESeq2_1.28.1 BiocFileCache_1.12.0
#> [5] dbplyr_1.4.4 GO.db_3.11.4
#> [7] org.Hs.eg.db_3.11.4 AnnotationDbi_1.50.2
#> [9] airway_1.8.0 SummarizedExperiment_1.18.2
#> [11] DelayedArray_0.14.1 matrixStats_0.56.0
#> [13] Biobase_2.48.0 GenomicRanges_1.40.0
#> [15] GenomeInfoDb_1.24.2 IRanges_2.22.2
#> [17] S4Vectors_0.26.1 BiocGenerics_0.34.0
#> [19] BiocSet_1.2.2 dplyr_1.0.0
#> [21] BiocStyle_2.16.0
#>
#> loaded via a namespace (and not attached):
#> [1] httr_1.4.2 splines_4.0.2 bit64_0.9-7.1
#> [4] assertthat_0.2.1 BiocManager_1.30.10 blob_1.2.1
#> [7] GenomeInfoDbData_1.2.3 Rsamtools_2.4.0 yaml_2.2.1
#> [10] pillar_1.4.6 RSQLite_2.2.0 lattice_0.20-41
#> [13] glue_1.4.1 digest_0.6.25 RColorBrewer_1.1-2
#> [16] XVector_0.28.0 colorspace_1.4-1 htmltools_0.5.0
#> [19] Matrix_1.2-18 plyr_1.8.6 XML_3.99-0.4
#> [22] pkgconfig_2.0.3 genefilter_1.70.0 bookdown_0.20
#> [25] zlibbioc_1.34.0 xtable_1.8-4 purrr_0.3.4
#> [28] scales_1.1.1 BiocParallel_1.22.0 annotate_1.66.0
#> [31] KEGGREST_1.28.0 ggplot2_3.3.2 generics_0.0.2
#> [34] ellipsis_0.3.1 cli_2.0.2 survival_3.2-3
#> [37] magrittr_1.5 crayon_1.3.4 memoise_1.1.0
#> [40] evaluate_0.14 fansi_0.4.1 tools_4.0.2
#> [43] formatR_1.7 lifecycle_0.2.0 stringr_1.4.0
#> [46] locfit_1.5-9.4 munsell_0.5.0 Biostrings_2.56.0
#> [49] compiler_4.0.2 rlang_0.4.7 grid_4.0.2
#> [52] RCurl_1.98-1.2 rappdirs_0.3.1 bitops_1.0-6
#> [55] rmarkdown_2.3 gtable_0.3.0 DBI_1.1.0
#> [58] curl_4.3 R6_2.4.1 GenomicAlignments_1.24.0
#> [61] knitr_1.29 rtracklayer_1.48.0 bit_1.1-15.2
#> [64] utf8_1.1.4 stringi_1.4.6 Rcpp_1.0.5
#> [67] geneplotter_1.66.0 vctrs_0.3.2 png_0.1-7
#> [70] tidyselect_1.1.0 xfun_0.15