rTRM: an R package for the identification of transcription regulatory modules (TRMs)
1 May 2024
Package
rTRM 1.43.0
Contents
1 Introduction
Transcription factors (TFs) bind to short motifs in the DNA and regulate the expression of target genes in a cell type and time dependent fashion. TFs do so by cooperating with other TFs in what it is called Transcriptional Regulatory Modules (TRMs). These TRMs contain different TFs and form a combinatorial code that explains TF specificity. We have implemented a method for the identification of TRMs that integrates information about binding locations from a single ChIP-seq experiment, computational estimation of TF binding, gene expression and protein-protein interaction (PPI) data (Diez, Hutchins, and Miranda-Saavedra 2014); see workflow figure). rTRM implements the methods required for the integration of PPI information (step 4 in workflow). To do so, rTRM tries to identify TFs that are bound to a target TF (the one with experimental evidence- i.e. ChIP-seq data) either directly or through a bridge protein. This package has been used to identify cell-type independent and dependent TRMs associated with Stat3 functions (Hutchins et al. 2013). Also, it has been used to identify TRMs in embryonic and hematopoietic stem cells as part of the publication presenting the methodology (Diez, Hutchins, and Miranda-Saavedra 2014). Here we present the basic capabilities of rTRM with a naive example, a case study showing the identification of Sox2 related TRM in ESCs as performed in the paper describing rTRM (Diez, Hutchins, and Miranda-Saavedra 2014), and a complete workflow in R using the PWMEnrich package for the motif enrichment step.
Figure 1: Workflow for the identification of TRMs
Steps 1-3 can be performed with standard Bioconductor approaches. rTRM implements a method to perform step 4.
2 Minimal example
In this minimal example a dummy network is search to identify TRMs focused around a target node, N6, with query nodes being N7, N12 and N28. By default findTRM find nodes separated a max distance of 0 (i.e. nodes directly connected). We change this with parameter max.bridge = 1. Because node N28 is separated by two other nodes from the target node N6, it is not included in the predicted TRM. By default findTRM returns an object of class igraph, which can be used with plotTRM, plotTRMlegend and other rTRM functions. However, it is possible to directly obtain a graphNEL object (from the Bioconductor package graph), setting type = graphNEL. Of course it is possible to also use the igraph.to.graphNEL function in the igraph package to transform an igraph object into a graphNEL object.
# load the rTRM package
library(rTRM)
# load network example.
load(system.file(package = "rTRM", "extra/example.rda"))
# plot network
plot(g, vertex.size = 20, vertex.label.cex = .8, layout = layout.graphopt)## This graph was created by an old(er) igraph version.
## Call upgrade_graph() on it to use with the current igraph version
## For now we convert it on the fly...
Figure 2: Identification of a TRM 1from a test network (left)
In the resulting TRM (right) dark blue indicates the target node light blue are query nodes and white nodes are bridge nodes
# define target and query nodes:
target <- "N6"
query <- c("N7", "N12", "N28")
# find TRM:
s <- findTRM(g, target = target, query = query, method = "nsa", max.bridge = 1)## removing: 4 nodes out of 8 [keeping 4 nodes]# annotate nodes:
V(s)$color <- "white"
V(s)[ name %in% query]$color <- "steelblue2"
V(s)[ name %in% target]$color <- "steelblue4"
# plot:
plot(s,vertex.size=20,vertex.label.cex=.8)
Figure 3: Identification of a TRM 1from a test network (left)
In the resulting TRM (right) dark blue indicates the target node light blue are query nodes and white nodes are bridge nodes
3 Introduction to the rTRM package
rTRM relies on a series of optimizations. For example in the publication we used PWMs for vertebrate species compiled from different sources. This assumes the binding specificities of TFs will be conserved on all these species. Recent comparison between mouse and human PWMs suggests that this is the case for most TFs Jolma et al. (2013). rTRM also relies on protein-protein interaction data, and so provides utilities to download data from the BioGRID database (see below). As some of these functionalities are further integrated with existing Bioconductor functionality they may be defunct in the future.
3.1 Database
Information about TFs, including Position Specific Weight (PWMs) matrices, mapping to Entrez Gene identifiers, orthologs in mouse and human and other annotations are stored as a SQLite database. rTRM provides a basic API for accessing the data. Below there are some examples.
To obtain PWMs:
pwm <- getMatrices()
head(pwm, 1)## $MA0009.1
## 1 2 3 4 5 6 7 8 9 10 11
## A 0.05 0.025 1 0 0 0 0 0.00 0.025 1 0.775
## C 0.70 0.025 0 0 0 0 0 0.05 0.175 0 0.125
## G 0.20 0.000 0 1 1 0 1 0.00 0.700 0 0.000
## T 0.05 0.950 0 0 0 1 0 0.95 0.100 0 0.100To get annotations:
ann <- getAnnotations()
head(ann)## row_names pwm_id symbol family domain binding source
## 1 1 MA0009.1 T T monomer jaspar
## 2 2 MA0059.1 MYC::MAX Helix-Loop-Helix dimer jaspar
## 3 3 MA0146.1 Zfx BetaBetaAlpha-zinc finger monomer jaspar
## 4 4 MA0132.1 Pdx1 Homeo monomer jaspar
## 5 5 MA0162.1 Egr1 BetaBetaAlpha-zinc finger monomer jaspar
## 6 6 MA0093.1 USF1 Helix-Loop-Helix monomer jaspar
## note
## 1 2010
## 2 2010
## 3 2010
## 4 2010
## 5 2010
## 6 2010To get map of TFs to genes:
map <- getMaps()
head(map)## row_names pwm_id entrezgene organism
## 1 1 MA0009.1 20997 mouse
## 2 2 MA0059.1 4609 human
## 3 3 MA0059.1 4149 human
## 4 4 MA0146.1 22764 mouse
## 5 5 MA0132.1 18609 mouse
## 6 6 MA0162.1 13653 mouseTo get map of TFs to ortholog genes:
o <- getOrthologs(organism = "mouse")
head(o)## row_names entrezgene organism map_entrezgene map_organism
## 1 775 20997 mouse 20997 mouse
## 2 776 4609 human 107771 mouse
## 3 777 4149 human 17187 mouse
## 4 778 22764 mouse 22764 mouse
## 5 779 18609 mouse 18609 mouse
## 6 780 13653 mouse 13653 mouseIt is possible to map motif ids to entrezgene ids in the target organism (only between human and mouse). This is useful when all the information about existing PWMs is desired, as some TF binding affinities have only been studied in one organism.
getOrthologFromMatrix("MA0009.1", organism="human")## [1] "6862"getOrthologFromMatrix("MA0009.1", organism="mouse")## [1] "20997"3.2 Interactome data
rTRM requires information about protein-protein interactions (PPIs) for its predictions and includes interactome (PPI network) data from the BioGRID database (???). Currently mouse and human interactomes are supported. The networks are provided as an igraph object. To access the data use:
# check statistics about the network.
biogrid_mm()## Mouse PPI network data of class igraph## This graph was created by an old(er) igraph version.
## Call upgrade_graph() on it to use with the current igraph version
## For now we convert it on the fly...## Number of nodes: 5315## Number of edges: 11695## Source: The BioGRID (http://www.thebiogrid.org)## Release: 3.4.128## Downloaded: 2015-09-17## Use data(biogrid_mm) to load it# load mouse PPI network:
data(biogrid_mm)The amount of available PPI data increases rapidly so it is desirable to have a way to access the newest data conveniently. rTRM includes support for direct download and processing of PPI data from the BioGRID database. The PPI network is stored as an igraph object that can be readily used with rTRM or stored for later use. Below there is an example of the BioGRID database update procedure.
# obtain dataset.
db <- getBiogridData() # retrieves latest release.
# db = getBiogridData("3.2.96") # to get a specific release.
# check release:
db$release
db$data
# process PPI data for different organisms (currently supported human and mouse):
biogrid_hs <- processBiogrid(db, org = "human")
biogrid_mm <- processBiogrid(db, org = "mouse")PPI data from other databases could be used as long as it is formatted as an igraph object with the name attribute containing entrezgene identifiers and the label attribute containing the symbol.
3.3 Using PSICQUIC package to obtain protein-protein interactions
One possibility available from Bioconductor is to use the package PSICQUIC to obtain PPI data. PSICQUIC provides access to different databases of PPIs, including BioGRID and STRINGS, and databases of cellular networks like KEGG or Reactome. For example, to obtain the human BioGRID data (NOTE: named BioGrid in PSICQUIC):
library(PSICQUIC)
psicquic <- PSICQUIC()
providers(psicquic)
# obtain BioGrid human PPIs (as data.frame):
tbl <- interactions(psicquic, species="9606",provider="BioGrid")
# the target and source node information needs to be polished (i.e. must be Entrez gene id only)
biogrid_hs <- data.frame(source=tbl$A,target=tbl$B)
biogrid_hs$source <- sub(".*locuslink:(.*)\\|BIOGRID:.*","\\1", biogrid_hs$source)
biogrid_hs$target <- sub(".*locuslink:(.*)\\|BIOGRID:.*","\\1", biogrid_hs$target)
# create graph.
library(igraph)
biogrid_hs <- graph.data.frame(biogrid_hs,directed=FALSE)
biogrid_hs <- simplify(biogrid_hs)
# annotate with symbols.
library(org.Hs.eg.db)
V(biogrid_hs)$label <- select(org.Hs.eg.db,keys=V(biogrid_hs)$name,columns=c("SYMBOL"))$SYMBOL4 Case study: TRM associated with Sox2 in embryonic stem cells (ESCs)
Sox2 is a TF involved in the determination and maintainance of pluripotency in embryonic stem cells (ESCs). Sox2 forms a transcriptional regulatory module with Nanog and Pou5f1 (Oct4), and together determine ESCs phenotype. Other TFs important to this process are Erssb and Klf4. In this case study we want to identify TRMs associated with Sox2. ChIP-seq data for Sox2 was obtained from Chen et al. (2008) and motif enrichment analysis performed with HOMER Heinz et al. (2010), followed by matching against our library of PWMs using TOMTOM Gupta et al. (2007). The starting dataset is the TOMTOM output file with the motifs enriched in the Sox2 binding regions.
# read motif enrichment results.
motif_file <- system.file("extra/sox2_motif_list.rda", package = "rTRM")
load(motif_file)
length(motif_list)## [1] 177head(motif_list)## [1] "MA0039.2" "MA0071.1" "MA0075.1" "MA0077.1" "MA0078.1" "MA0112.2"First, we read the motifs and convert them into gene identifiers (i.e. Entrez Gene identifier). To do this we use the function getOrthologFromMatrix, which takes a list of motif identifiers and the target organism as parameters. The function returns a list with the Entrez Gene ids.
# get the corresponding gene.
tfs_list <- getOrthologFromMatrix(motif_list, organism = "mouse")
tfs_list <- unique(unlist(tfs_list, use.names = FALSE))
length(tfs_list)## [1] 98head(tfs_list)## [1] "18609" "20682" "13983" "20665" "18291" "18227"Next, we need a list of genes expressed in ESC. For this, the dataset was obtained from GEO (GSE27708; Ho et al. (2011)) and processed using the custom CDFs from the BrainArray project Dai et al. (2005) and the rma function from the package affy Gautier et al. (2004). Genes not expressed were filtered by removing all genes with log2 expression < 5 in all samples.
# load expression data.
eg_esc_file <- system.file("extra/ESC-expressed.txt", package = "rTRM")
eg_esc <- scan(eg_esc_file, what = "")
length(eg_esc)## [1] 8734head(eg_esc)## [1] "100008567" "100017" "100019" "100037258" "100038489" "100039781"tfs_list_esc <- tfs_list[tfs_list %in% eg_esc]
length(tfs_list_esc)## [1] 22head(tfs_list_esc)## [1] "26380" "18999" "20674" "16600" "26379" "13984"Next, we load the PPI network and filter out potential degree outliers and proteins not expressed in the paired expression data.
# load and process PPI data.
biogrid_mm()## Mouse PPI network data of class igraph## This graph was created by an old(er) igraph version.
## Call upgrade_graph() on it to use with the current igraph version
## For now we convert it on the fly...## Number of nodes: 5315## Number of edges: 11695## Source: The BioGRID (http://www.thebiogrid.org)## Release: 3.4.128## Downloaded: 2015-09-17## Use data(biogrid_mm) to load itdata(biogrid_mm)
ppi <- biogrid_mm
vcount(ppi)## This graph was created by an old(er) igraph version.
## Call upgrade_graph() on it to use with the current igraph version
## For now we convert it on the fly...## [1] 5315ecount(ppi)## [1] 11695# remove outliers.
f <- c("Ubc", "Sumo1", "Sumo2", "Sumo3")
f <- select(org.Mm.eg.db, keys = f, columns = "ENTREZID", keytype = "SYMBOL")$ENTREZID## 'select()' returned 1:1 mapping between keys and columnsf## [1] "22190" "22218" "170930" "20610"ppi <- removeVertices(ppi, f)
vcount(ppi)## [1] 4984ecount(ppi)## [1] 11081# filter by expression.
ppi_esc <- induced.subgraph(ppi, V(ppi)[ name %in% eg_esc ])## Warning: `induced.subgraph()` was deprecated in igraph 2.0.0.
## ℹ Please use `induced_subgraph()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.vcount(ppi_esc)## [1] 3109ecount(ppi_esc)## [1] 4889# ensure a single component.
ppi_esc <- getLargestComp(ppi_esc)
vcount(ppi_esc)## [1] 2576ecount(ppi_esc)## [1] 4851To identify TRMs we define a target TF (the one the ChIP-seq data comes from) and some query TFs (the ones with enriched binding sites in the neighborhood of the target TF).
# define target.
target <- select(org.Mm.eg.db,keys="Sox2",columns="ENTREZID",keytype="SYMBOL")$ENTREZID## 'select()' returned 1:1 mapping between keys and columnstarget## [1] "20674"# find TRM.
s <- findTRM(ppi_esc, target, tfs_list_esc, method = "nsa", max.bridge = 1)## 11 query nodes NOT FOUND in network-- removed## removing: 320 nodes out of 328 [keeping 8 nodes]vcount(s)## [1] 8ecount(s)## [1] 15Finally, we layout the network using a customized concentric layout and plot the network and the legend.
# generate layout (order by cluster, then label)
cl <- getConcentricList(s, target, tfs_list_esc)
l <- layout.concentric(s, cl, order = "label")
# plot TRM.
plotTRM(s, layout = l, vertex.cex = 15, label.cex = .8)