Contents

1 Introduction

xcore package provides a framework for transcription factor (TF) activity modeling based on their molecular signatures and user’s gene expression data. Additionally, xcoredata package provides a collection of pre-processed TF molecular signatures constructed from ChiP-seq experiments available in ReMap2020 or ChIP-Atlas databases.

xcore flow chart

xcore flow chart

xcore use ridge regression to model changes in expression as a linear combination of molecular signatures and find their unknown activities. Obtained, estimates can be further tested for significance to select molecular signatures with the highest predicted effect on the observed expression changes.

2 Installation

Currently, you can install xcore package for R 4.1 from github using:

devtools::install_github("bkaczkowski/xcore@r4_1")

xcore package is currently being submitted for mid April release of the Bioconductor. After mid April, it will be possible to install it as Bioconductor package.

BiocManager::install("xcore") # after mid April

3 Gene expression modeling in context of rinderpest infection

library("xcore")

Here we will use a subset of 293SLAM rinderpest infection dataset from FANTOM5. This subset contains expression counts for 0, 12 and 24 hours post infection samples and only for a subset of FANTOM5 promoters. We can find this example data shipped together with xcore package.

This is a simple count matrix, with columns corresponding to samples and rows corresponding to FANTOM5 promoters.

data("rinderpest_mini")
head(rinderpest_mini)
##                                              00hr_rep1 00hr_rep2 00hr_rep3
## hg19::chr1:10003372..10003465,-;hg_10258.1          52        46        57
## hg19::chr1:10003486..10003551,+;hg_541.1            39        42        27
## hg19::chr1:100111580..100111773,+;hg_4181.1          1         0         2
## hg19::chr1:100232177..100232198,-;hg_13495.1        15         9        26
## hg19::chr1:100315613..100315691,+;hg_4187.1         95       109       110
## hg19::chr1:100435545..100435597,+;hg_4201.1        141       129       101
##                                              12hr_rep1 12hr_rep2 12hr_rep3
## hg19::chr1:10003372..10003465,-;hg_10258.1          54        50        53
## hg19::chr1:10003486..10003551,+;hg_541.1            35        30        40
## hg19::chr1:100111580..100111773,+;hg_4181.1          1         2         3
## hg19::chr1:100232177..100232198,-;hg_13495.1        20        16        13
## hg19::chr1:100315613..100315691,+;hg_4187.1        112        94       103
## hg19::chr1:100435545..100435597,+;hg_4201.1        132       106       125
##                                              24hr_rep1 24hr_rep2 24hr_rep3
## hg19::chr1:10003372..10003465,-;hg_10258.1          11        12        12
## hg19::chr1:10003486..10003551,+;hg_541.1            22        34        50
## hg19::chr1:100111580..100111773,+;hg_4181.1          0         1         1
## hg19::chr1:100232177..100232198,-;hg_13495.1         7        13        10
## hg19::chr1:100315613..100315691,+;hg_4187.1         43        74        89
## hg19::chr1:100435545..100435597,+;hg_4201.1         84       100       121

First we need to construct a design matrix describing our experiment design.

design <- matrix(
  data = c(1, 0, 0,
           1, 0, 0,
           1, 0, 0,
           0, 1, 0,
           0, 1, 0,
           0, 1, 0,
           0, 0, 1,
           0, 0, 1,
           0, 0, 1),
  ncol = 3,
  nrow = 9,
  byrow = TRUE,
  dimnames = list(
    c(
      "00hr_rep1",
      "00hr_rep2",
      "00hr_rep3",
      "12hr_rep1",
      "12hr_rep2",
      "12hr_rep3",
      "24hr_rep1",
      "24hr_rep2",
      "24hr_rep3"
    ),
    c("00hr", "12hr", "24hr")
  )
)

print(design)
##           00hr 12hr 24hr
## 00hr_rep1    1    0    0
## 00hr_rep2    1    0    0
## 00hr_rep3    1    0    0
## 12hr_rep1    0    1    0
## 12hr_rep2    0    1    0
## 12hr_rep3    0    1    0
## 24hr_rep1    0    0    1
## 24hr_rep2    0    0    1
## 24hr_rep3    0    0    1

Next, we need to preprocess the counts using prepareCountsForRegression function. Here, CAGE expression tags for each sample are filtered for lowly expressed promoters, normalized for the library size and transformed into counts per million (CPM). Finally, CPM are log2 transformed with addition of pseudo count 1. Moreover, we designate the base level samples, 0 hours after treatment in our example, from which basal expression level is calculated. This basal level will be used as a reference when modeling the expression changes.

mae <- prepareCountsForRegression(counts = rinderpest_mini,
                                  design = design,
                                  base_lvl = "00hr")

xcore models the expression as a function of molecular signatures. Such signatures can be constructed eg. from the known transcription factor binding sites (see Constructing molecular signatures section). Here, we will take advantage of pre-computed molecular signatures found in the xcoredata package. Particularly, molecular signatures constructed from ReMap2020 against FANTOM5 annotation.

Molecular signatures can be conveniently accessed using the ExperimentHub interface. This functionality will be available after mid April.

library("ExperimentHub") # after mid April

eh <- ExperimentHub()
query(eh, "xcoredata")
## ExperimentHub with 6 records
## # snapshotDate(): 2022-04-19
## # $dataprovider: ReMap, RIKEN, ChIP-Atlas
## # $species: Homo sapiens, NA
## # $rdataclass: dgCMatrix, data.table, GRanges
## # additional mcols(): taxonomyid, genome, description,
## #   coordinate_1_based, maintainer, rdatadateadded, preparerclass, tags,
## #   rdatapath, sourceurl, sourcetype 
## # retrieve records with, e.g., 'object[["EH7298"]]' 
## 
##            title                  
##   EH7298 | chip_atlas_meta        
##   EH7299 | remap_meta             
##   EH7300 | chip_atlas_promoters_f5
##   EH7301 | remap_promoters_f5     
##   EH7302 | promoters_f5           
##   EH7303 | promoters_f5_core
remap_promoters_f5 <- eh[["EH7301"]] # after mid April

Currently, molecular signatures can be accessed directly from xcoredata package.

devtools::install_gitlab("mcjmigdal/xcoredata@r4_1")
data("remap_promoters_f5", package = "xcoredata")

Molecular signature is a simple binary matrix indicating if a transcription factor binding site was found or not found in promoter vicinity.

print(remap_promoters_f5[3:6, 3:6])
## 4 x 4 sparse Matrix of class "dgCMatrix"
##                                             MLLT1.GM12878.ENCSR552XSN
## hg19::chr1:10003372..10003465,-;hg_10258.1                          1
## hg19::chr1:10003486..10003551,+;hg_541.1                            1
## hg19::chr1:100110807..100110818,+;hg_4172.1                         1
## hg19::chr1:100110851..100110860,+;hg_4173.1                         1
##                                             ZBTB10.HEK293.ENCSR004PLU
## hg19::chr1:10003372..10003465,-;hg_10258.1                          1
## hg19::chr1:10003486..10003551,+;hg_541.1                            1
## hg19::chr1:100110807..100110818,+;hg_4172.1                         .
## hg19::chr1:100110851..100110860,+;hg_4173.1                         .
##                                             ZFP3.HEK293.ENCSR134QIE
## hg19::chr1:10003372..10003465,-;hg_10258.1                        .
## hg19::chr1:10003486..10003551,+;hg_541.1                          .
## hg19::chr1:100110807..100110818,+;hg_4172.1                       .
## hg19::chr1:100110851..100110860,+;hg_4173.1                       .
##                                             ZNF2.HEK293.ENCSR011CKE
## hg19::chr1:10003372..10003465,-;hg_10258.1                        1
## hg19::chr1:10003486..10003551,+;hg_541.1                          1
## hg19::chr1:100110807..100110818,+;hg_4172.1                       .
## hg19::chr1:100110851..100110860,+;hg_4173.1                       .

3.0.1 Computational resources consideration

Running xcore using extensive ReMap2020 or ChIP-Atlas molecular signatures, can be quite time and memory consuming. As an example, modeling the example 293SLAM rinderpest infection dataset with ReMap2020 signatures matrix took around 10 minutes to compute and used up to 8 GB RAM on Intel(R) Xeon(R) CPU E5-2680 v3 using 2 cores. This unfortunately exceeds the resources available for Bioconductor vignettes. This being said, we will further proceed by taking only a subset of ReMap2020 signatures such that we fit into the time limits. However, for running xcore in a normal setting the whole molecular signatures matrix should be used!

# here we subset ReMap2020 molecular signatures matrix for the purpose of the
# vignette only! In a normal setting the whole molecular signatures matrix should
# be used!
set.seed(432123)
i <- sample(x = seq_len(ncol(remap_promoters_f5)), size = 100, replace = FALSE)
remap_promoters_f5 <- remap_promoters_f5[, i]

To add signatures to our MultiAssayExperiment object we can use addSignatures function. As you add your signatures remember to give them unique names.

mae <- addSignatures(mae, remap = remap_promoters_f5)

When we examine newly added signatures, we can see that some of them does not overlap any of the promoters. On the other side, depending on the signatures quality, we could expect to see signatures that overlap all of the promoters.

summary(colSums(mae[["remap"]]))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     1.0    99.5  1300.5  2499.3  4148.0  9561.0

While, we do not provide detailed guidelines to signatures filtering, filtering out at least signatures overlapping no or all promoters is mandatory. Here, we filter signatures that overlap less than 5% or more than 95% of promoters using filterSignatures function.

mae <- filterSignatures(mae, min = 0.05, max = 0.95)

At this stage we can investigate mae to see how many signatures are available after intersecting with counts and filtering. number of columns in signatures corresponds to number of molecular signatures which will be used to build a model.

print(mae)
## A MultiAssayExperiment object of 3 listed
##  experiments with user-defined names and respective classes.
##  Containing an ExperimentList class object of length 3:
##  [1] U: matrix with 10388 rows and 1 columns
##  [2] Y: matrix with 10388 rows and 6 columns
##  [3] remap: dgCMatrix with 10388 rows and 63 columns
## Functionality:
##  experiments() - obtain the ExperimentList instance
##  colData() - the primary/phenotype DataFrame
##  sampleMap() - the sample coordination DataFrame
##  `$`, `[`, `[[` - extract colData columns, subset, or experiment
##  *Format() - convert into a long or wide DataFrame
##  assays() - convert ExperimentList to a SimpleList of matrices
##  exportClass() - save data to flat files

Finally, we can run our expression modeling using modelGeneExpression function. modelGeneExpression can take advantage of parallelization to speed up the calculations. To use it we need to register a parallel backend.

3.0.2 Parallel computing

While there are many parallel backends to choose from, internally xcore uses foreach to implement parallel computing. Having this in mind we should use a backend supported by foreach.

Here we are using doParallel backend, together with BiocParallel package providing unified interface across different OS. Those packages can be installed with:

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("BiocParallel")
install.packages("doParallel")

Our modeling will inherently contain some level of randomness due to using cross validation, here we set the seed so that our results can be replicated.

This step is time consuming as we need to train a separate model for each of our samples. To speed up the calculations we will lower number of folds used in cross-validation procedure from 10 to 5 using nfolds argument, which is internally passed to cv.glmnet function. Moreover here we only use a subset of the ReMap2020 molecular signatures matrix.

# register parallel backend
doParallel::registerDoParallel(2L)
BiocParallel::register(BiocParallel::DoparParam(), default = TRUE)

# set seed
set.seed(314159265)

res <- modelGeneExpression(
  mae = mae,
  xnames = "remap",
  nfolds = 5)

Output returned by modelGeneExpression is a list with following elements:

  • regression_models - a list holding cv.glmnet objects for each sample.
  • pvalues - list of data.frames for each sample, each holding signatures estimates, their estimated standard errors and p-values.
  • zscore_avg - a matrix holding replicate average molecular signatures Z-scores with columns corresponding to groups in the design.
  • coef_avg - a matrix holding replicate average molecular signatures activities with columns corresponding to groups in the design.
  • results - a data.frame holding replicate average molecular signatures and overall molecular signatures Z-score calculated over all groups in the design.

results is likely of most interest. As you can see below this data.frame holds replicate averaged molecular signatures activities, as well as overall Z-score which can be used to rank signatures.

head(res$results$remap)
##                                 name        12hr        24hr  z_score
## 17            MYC.NCI-H2171.GSE36354 -0.12198729 -0.20038927 30.07425
## 51      E2F4.lymphoblastoid.GSE21488  0.04584327  0.20138692 23.27914
## 10            ATF3.K-562.ENCSR632DCH  0.03166150  0.18521027 19.93733
## 28 NELFA.HeLa_Flavo-0-H2O2.GSE100742 -0.03904617 -0.14593124 13.83496
## 31           MXD4.Hep-G2.ENCSR441KFW -0.07144238 -0.07360724 10.62700
## 59              MYC.HFF_OHT.GSE65544  0.03300246 -0.06842935 10.24881

To visualize the result we can plot a heatmap of molecular signatures activities for the top identified molecular signatures.

top_signatures <- head(res$results$remap, 10)
library(pheatmap)
pheatmap::pheatmap(
  mat = top_signatures[, c("12hr", "24hr")],
  scale = "none",
  labels_row = top_signatures$name,
  cluster_cols = FALSE,
  color = colorRampPalette(c("blue", "white", "red"))(15),
  breaks = seq(from = -0.1, to = 0.1, length.out = 16),
  main = "SLAM293 molecular signatures activities"
)
Predicted activities for the top-scoring molecular signatures

Figure 1: Predicted activities for the top-scoring molecular signatures

As we only used a random subset of ReMap2020 molecular signatures the obtained result most likely has no biological meaning. Nonetheless, some general comment on results interpretation can be made. The top-scoring transcription factors can be hypothesized to be associated with observed changes in gene expression. Moreover, the predicted activities indicate if promoters targeted by a specific signature tend to be downregulated or upregulated.

4 Constructing molecular signatures

Constructing molecular signatures is as simple as intersecting promoters annotation with transcription factors binding sites (TFBS). For example, here we use FANTOM5’ promoters annotation found in the xcoredata package and predicted CTCF binding sites available in CTCF AnnotationHub’s resource.

Promoter annotations and TFBS should be stored as a GRanges object. Moreover, it is required that the promoter/TF name is held under the name column. Next we can construct a molecular signatures matrix using getInteractionMatrix function, where we can specify the ext parameter that will control how much the promoters are extended in both directions before intersecting with TFBS.

# obtain promoters annotation; *promoter name must be stored under column 'name'
promoters_f5 <- eh[["EH7303"]]
head(promoters_f5)
## GRanges object with 6 ranges and 5 metadata columns:
##       seqnames            ranges strand |                   name     score
##          <Rle>         <IRanges>  <Rle> |            <character> <numeric>
##   [1]     chr1   9943315-9943407      - | hg19::chr1:10003372...    102331
##   [2]     chr1   9943429-9943493      + | hg19::chr1:10003486...     69018
##   [3]     chr1 99646025-99646217      + | hg19::chr1:100111580..     87706
##   [4]     chr1 99766622-99766642      - | hg19::chr1:100232177..     25838
##   [5]     chr1 99850058-99850135      + | hg19::chr1:100315613..     92180
##   [6]     chr1 99969990-99970041      + | hg19::chr1:100435545..    180010
##       gene_type_gencode    ENTREZID      SYMBOL
##                <factor> <character> <character>
##   [1]    protein_coding       84328        LZIC
##   [2]    protein_coding       64802      NMNAT1
##   [3]    protein_coding       54873       PALMD
##   [4]    protein_coding      391059       FRRS1
##   [5]    protein_coding         178         AGL
##   [6]    protein_coding       23443     SLC35A3
##   -------
##   seqinfo: 25 sequences (1 circular) from hg38 genome
# obtain predicted CTCF binding for hg38;
# *TF/motif name must be stored under column 'name'
library("AnnotationHub")
ah <- AnnotationHub()
CTCF_hg38 <- ah[["AH95566"]]
CTCF_hg38$name <- CTCF_hg38$motif
head(CTCF_hg38)
## GRanges object with 6 ranges and 6 metadata columns:
##       seqnames        ranges strand |       motif     score   p.value   q.value
##          <Rle>     <IRanges>  <Rle> | <character> <numeric> <numeric> <numeric>
##   [1]     chr1   11223-11241      - |    MA0139.1   24.4754  1.34e-09    0.0216
##   [2]     chr1   11281-11299      - |    MA0139.1   22.7377  1.01e-08    0.0398
##   [3]     chr1   24782-24800      - |    MA0139.1   17.3770  7.11e-07    0.2350
##   [4]     chr1   91420-91438      + |    MA0139.1   16.2951  1.41e-06    0.3080
##   [5]     chr1 104985-105003      - |    MA0139.1   16.7869  1.04e-06    0.2750
##   [6]     chr1 132815-132833      - |    MA0139.1   15.3770  2.44e-06    0.3510
##                  sequence        name
##               <character> <character>
##   [1] TCGCCAGCAGGGGGCGCCC    MA0139.1
##   [2] GCGCCAGCAGGGGGCGCTG    MA0139.1
##   [3] CGTCCAGCAGATGGCGGAT    MA0139.1
##   [4] GTGGCACCAGGTGGCAGCA    MA0139.1
##   [5] CCAACAGCAGGTGGCAGCC    MA0139.1
##   [6] GCACCAGAAGGTGGCAGGA    MA0139.1
##   -------
##   seqinfo: 24 sequences from hg38 genome
# construct molecular signatures matrix
molecular_signature <- getInteractionMatrix(
  a = promoters_f5,
  b = CTCF_hg38,
  ext = 500
)
head(molecular_signature)
## 6 x 1 sparse Matrix of class "dgCMatrix"
##                                              MA0139.1
## hg19::chr1:10003372..10003465,-;hg_10258.1          .
## hg19::chr1:10003486..10003551,+;hg_541.1            1
## hg19::chr1:100111580..100111773,+;hg_4181.1         .
## hg19::chr1:100232177..100232198,-;hg_13495.1        .
## hg19::chr1:100315613..100315691,+;hg_4187.1         .
## hg19::chr1:100435545..100435597,+;hg_4201.1         .
sessionInfo()
## R version 4.2.0 RC (2022-04-19 r82224)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.4 LTS
## 
## Matrix products: default
## BLAS:   /home/biocbuild/bbs-3.15-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.15-bioc/R/lib/libRlapack.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_GB              LC_COLLATE=C              
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] GenomicRanges_1.48.0 GenomeInfoDb_1.32.0  IRanges_2.30.0      
##  [4] S4Vectors_0.34.0     pheatmap_1.0.12      glmnet_4.1-4        
##  [7] Matrix_1.4-1         xcoredata_0.99.2     ExperimentHub_2.4.0 
## [10] AnnotationHub_3.4.0  BiocFileCache_2.4.0  dbplyr_2.1.1        
## [13] BiocGenerics_0.42.0  xcore_1.0.0          BiocStyle_2.24.0    
## 
## loaded via a namespace (and not attached):
##  [1] bitops_1.0-7                  matrixStats_0.62.0           
##  [3] bit64_4.0.5                   RColorBrewer_1.1-3           
##  [5] doParallel_1.0.17             filelock_1.0.2               
##  [7] httr_1.4.2                    tools_4.2.0                  
##  [9] bslib_0.3.1                   utf8_1.2.2                   
## [11] R6_2.5.1                      colorspace_2.0-3             
## [13] DBI_1.1.2                     withr_2.5.0                  
## [15] tidyselect_1.1.2              bit_4.0.4                    
## [17] curl_4.3.2                    compiler_4.2.0               
## [19] cli_3.3.0                     Biobase_2.56.0               
## [21] DelayedArray_0.22.0           bookdown_0.26                
## [23] sass_0.4.1                    scales_1.2.0                 
## [25] rappdirs_0.3.3                stringr_1.4.0                
## [27] digest_0.6.29                 rmarkdown_2.14               
## [29] XVector_0.36.0                pkgconfig_2.0.3              
## [31] htmltools_0.5.2               MatrixGenerics_1.8.0         
## [33] highr_0.9                     fastmap_1.1.0                
## [35] limma_3.52.0                  rlang_1.0.2                  
## [37] RSQLite_2.2.12                shiny_1.7.1                  
## [39] shape_1.4.6                   jquerylib_0.1.4              
## [41] generics_0.1.2                jsonlite_1.8.0               
## [43] BiocParallel_1.30.0           dplyr_1.0.8                  
## [45] RCurl_1.98-1.6                magrittr_2.0.3               
## [47] GenomeInfoDbData_1.2.8        munsell_0.5.0                
## [49] Rcpp_1.0.8.3                  fansi_1.0.3                  
## [51] lifecycle_1.0.1               stringi_1.7.6                
## [53] yaml_2.3.5                    edgeR_3.38.0                 
## [55] SummarizedExperiment_1.26.0   zlibbioc_1.42.0              
## [57] grid_4.2.0                    blob_1.2.3                   
## [59] parallel_4.2.0                promises_1.2.0.1             
## [61] crayon_1.5.1                  lattice_0.20-45              
## [63] Biostrings_2.64.0             splines_4.2.0                
## [65] KEGGREST_1.36.0               magick_2.7.3                 
## [67] locfit_1.5-9.5                knitr_1.38                   
## [69] pillar_1.7.0                  codetools_0.2-18             
## [71] glue_1.6.2                    BiocVersion_3.15.2           
## [73] evaluate_0.15                 BiocManager_1.30.17          
## [75] MultiAssayExperiment_1.22.0   vctrs_0.4.1                  
## [77] png_0.1-7                     httpuv_1.6.5                 
## [79] foreach_1.5.2                 gtable_0.3.0                 
## [81] purrr_0.3.4                   assertthat_0.2.1             
## [83] cachem_1.0.6                  xfun_0.30                    
## [85] mime_0.12                     xtable_1.8-4                 
## [87] later_1.3.0                   survival_3.3-1               
## [89] tibble_3.1.6                  iterators_1.0.14             
## [91] AnnotationDbi_1.58.0          memoise_2.0.1                
## [93] ellipsis_0.3.2                interactiveDisplayBase_1.34.0