4 ArchR workflow and different weight methods
In this chapter, we illustrate the epiregulon workflow starting from an ArchR project and compare the different weight estimation methods.
The dataset consists of unpaired scATACseq/scRNAseq of parental LNCaP cells treated with DMSO, Enzalutamide and Enza resistant cells. The dataset was taken from Taavitsainen et al GSE168667 and GSE168668.
4.1 Data preparation
Please refer to the full ArchR manual for instructions
Before running Epiregulon, the following analyses need to be completed:
- Obtain a peak matrix on scATACseq by using addGroupCoverages > addReproduciblePeakSet > addPeakMatrix. See chapter 10 from ArchR manual
- RNA-seq integration.
- Perform dimensionality reduction from with either single modalities or joint scRNAseq and scATACseq using addCombinedDims
4.2 Load ArchR project
library(ArchR)
archR_project_path <- "/gstore/project/lineage/prostate/GSE168667/OUTPUT/multiome/"
proj <- loadArchRProject(path = archR_project_path, showLogo = FALSE)We verify that “GeneExpressionMatrix” and “PeakMatrix” are present for this tutorial.
## [1] "GeneIntegrationMatrix" "GeneScoreMatrix" "MotifMatrix" "PeakMatrix" "TileMatrix"We will use the joint reducedDims - “LSI_Combined” and joint embeddings - “UMAP_Combined”
## LSI1 LSI2 LSI3 LSI4 LSI5
## SRR13927735#TTATGTCTCCAGGTAT-1 -2.713935 -0.3677949 -0.4484238 -0.30645138 -0.046845365
## SRR13927735#TATTGCTCATCAGAAA-1 -2.642781 -0.2767556 -0.9142714 -0.19675812 0.075746940
## SRR13927735#TTCGATTGTAGGGTTG-1 -2.322865 -0.1543080 -1.4106049 -0.08891276 0.019873276
## SRR13927735#CATTCATTCGGATGTT-1 -2.572976 -0.1917188 -1.0464294 -0.12660121 0.009947438
## SRR13927735#ACGTTAGGTCAACTGT-1 -2.478552 -0.1776639 -1.1037295 -0.22976613 -0.150097539
## SRR13927735#AAATGCCCAGCAATGG-1 -2.595352 -0.3803464 -0.7770309 -0.52431765 -0.243074591## iLSI_Combined#UMAP_Dimension_1 iLSI_Combined#UMAP_Dimension_2
## SRR13927735#TTATGTCTCCAGGTAT-1 -9.622903 -0.2908237
## SRR13927735#TATTGCTCATCAGAAA-1 -9.360211 -0.2892935
## SRR13927735#TTCGATTGTAGGGTTG-1 -8.617347 -0.2154103
## SRR13927735#CATTCATTCGGATGTT-1 -9.285448 -0.3267481
## SRR13927735#ACGTTAGGTCAACTGT-1 -8.809260 -0.2168703
## SRR13927735#AAATGCCCAGCAATGG-1 -9.261216 0.32003564.3 Retrieve matrices from ArchR project
Retrieve gene expression and peak matrix from the ArchR project
GeneExpressionMatrix <- getMatrixFromProject(
ArchRProj = proj,
useMatrix = "GeneIntegrationMatrix",
useSeqnames = NULL,
verbose = TRUE,
binarize = FALSE,
threads = 1,
logFile = "x"
)## 2024-03-19 08:01:18.238282 : Organizing colData, 1.926 mins elapsed.## 2024-03-19 08:01:18.424896 : Organizing rowData, 1.93 mins elapsed.## 2024-03-19 08:01:18.429259 : Organizing rowRanges, 1.93 mins elapsed.## 2024-03-19 08:01:18.436044 : Organizing Assays (1 of 1), 1.93 mins elapsed.## 2024-03-19 08:01:28.057361 : Constructing SummarizedExperiment, 2.09 mins elapsed.## 2024-03-19 08:01:31.062914 : Finished Matrix Creation, 2.14 mins elapsed.PeakMatrix <- getMatrixFromProject(
ArchRProj = proj,
useMatrix = "PeakMatrix",
useSeqnames = NULL,
verbose = TRUE,
binarize = FALSE,
threads = 1,
logFile = "x"
)## 2024-03-19 08:03:00.665676 : Organizing colData, 1.486 mins elapsed.## 2024-03-19 08:03:01.440789 : Organizing rowData, 1.498 mins elapsed.## 2024-03-19 08:03:01.450369 : Organizing rowRanges, 1.499 mins elapsed.## 2024-03-19 08:03:01.462803 : Organizing Assays (1 of 1), 1.499 mins elapsed.## 2024-03-19 08:03:04.8355 : Constructing SummarizedExperiment, 1.555 mins elapsed.## 2024-03-19 08:03:27.161162 : Finished Matrix Creation, 1.927 mins elapsed.If we extract the gene expression from matrix, it will be in the form of RangedSummarizedExperiment. We can make use of ArchRMatrix2SCE to convert gene expression matrix to SingleCellExperiment object. It’s also important to note that gene expression from ArchR is library size normalized (not logged)
##
## Attaching package: 'epiregulon.archr'## The following objects are masked from 'package:epiregulon':
##
## addMotifScore, addTFMotifInfo, calculateP2G, getTFMotifInfoGeneExpressionMatrix <- ArchRMatrix2SCE(GeneExpressionMatrix)
rownames(GeneExpressionMatrix) <- rowData(GeneExpressionMatrix)$nameWe rename the assay name of the PeakMatrix as counts
Transfer embeddings from ArchR project to singleCellExperiment for visualization
reducedDim(GeneExpressionMatrix, "UMAP_Combined") <- getEmbedding(ArchRProj = proj,
embedding = "UMAP_Combined",
returnDF = TRUE)[colnames(GeneExpressionMatrix),]
# add cell label
GeneExpressionMatrix$label <- GeneExpressionMatrix$Cells
GeneExpressionMatrix$label[GeneExpressionMatrix$Treatment == "enzalutamide 48h"] <- "LNCaP–ENZ48"
GeneExpressionMatrix$label <- factor(GeneExpressionMatrix$label,
levels = c("LNCaP", "LNCaP–ENZ48", "LNCaP RES-A", "LNCaP RES-B"))Visualize singleCellExperiment by UMAP
scater::plotReducedDim(GeneExpressionMatrix,
dimred = "UMAP_Combined",
text_by = "label",
colour_by = "label")
4.4 Retrieve bulk TF ChIP-seq binding sites
First, we retrieve the information of TF binding sites collected from Cistrome and ENCODE ChIP-seq. Currently, human genomes HG19 and HG38 and mouse mm10 are available.
## snapshotDate(): 2024-03-11## see ?scMultiome and browseVignettes('scMultiome') for documentation## loading from cache## GRangesList object of length 1558:
## $AEBP2
## GRanges object with 2700 ranges and 0 metadata columns:
## seqnames ranges strand
## <Rle> <IRanges> <Rle>
## [1] chr1 9792-10446 *
## [2] chr1 942105-942400 *
## [3] chr1 984486-984781 *
## [4] chr1 3068932-3069282 *
## [5] chr1 3069411-3069950 *
## ... ... ... ...
## [2696] chrY 8465261-8465730 *
## [2697] chrY 11721744-11722260 *
## [2698] chrY 11747448-11747964 *
## [2699] chrY 19302661-19303134 *
## [2700] chrY 19985662-19985982 *
## -------
## seqinfo: 25 sequences from an unspecified genome; no seqlengths
##
## ...
## <1557 more elements>4.5 Link ATAC-seq peaks to target genes
Next, we compute peak to gene correlations using the addPeak2GeneLinks function from the ArchR package. The user would need
to supply a path to an ArchR project already containing peak and gene matrices, as well as Latent semantic indexing (LSI) dimensionality reduction.
# path to ArchR project
p2g <- calculateP2G(ArchR_path = archR_project_path,
useDim = "iLSI_Combined",
useMatrix = "GeneIntegrationMatrix",
threads = 1)## Setting ArchRLogging = FALSE## Using ArchR to compute peak to gene links...## 2024-03-19 08:03:55.593464 : Getting Available Matrices, 0 mins elapsed.## 2024-03-19 08:04:04.379045 : Filtered Low Prediction Score Cells (0 of 15522, 0), 0.015 mins elapsed.## 2024-03-19 08:04:05.042742 : Computing KNN, 0.026 mins elapsed.## 2024-03-19 08:04:06.92312 : Identifying Non-Overlapping KNN pairs, 0.058 mins elapsed.## 2024-03-19 08:04:09.163414 : Identified 498 Groupings!, 0.095 mins elapsed.## 2024-03-19 08:04:09.226792 : Getting Group RNA Matrix, 0.096 mins elapsed.## 2024-03-19 08:12:28.639651 : Getting Group ATAC Matrix, 8.419 mins elapsed.## 2024-03-19 08:17:42.039764 : Normalizing Group Matrices, 13.643 mins elapsed.## 2024-03-19 08:17:51.716899 : Finding Peak Gene Pairings, 13.804 mins elapsed.## 2024-03-19 08:17:52.303141 : Computing Correlations, 13.814 mins elapsed.## 2024-03-19 08:18:01.20975 : Completed Peak2Gene Correlations!, 13.962 mins elapsed.## DataFrame with 17979 rows and 8 columns
## idxATAC chr start end idxRNA target Correlation distance
## <integer> <factor> <integer> <integer> <integer> <character> <numeric> <numeric>
## 1 15 chr1 912762 913262 7 NOC2L 0.546722 46297
## 2 15 chr1 912762 913262 8 KLHL17 0.516539 47575
## 3 25 chr1 920261 920761 7 NOC2L 0.649425 38798
## 4 25 chr1 920261 920761 8 KLHL17 0.637711 40076
## 5 32 chr1 927728 928228 7 NOC2L 0.610240 31331
## ... ... ... ... ... ... ... ... ...
## 17975 210643 chrX 154542721 154543221 23496 CH17-340M24.3 0.611273 114492
## 17976 210643 chrX 154542721 154543221 23501 LAGE3 0.698033 63714
## 17977 210643 chrX 154542721 154543221 23506 IKBKG 0.518586 1716
## 17978 210643 chrX 154542721 154543221 23509 DKC1 0.526624 219771
## 17979 210665 chrX 154815200 154815700 23515 F8 0.547783 2114904.6 Add TF motif binding to peaks
The next step is to add the TF motif binding information by overlapping the regions of the peak matrix with the bulk chip-seq database. The user can supply an archR project path and this function will retrieve the peak matrix, or a peakMatrix in the form of a Granges object or RangedSummarizedExperiment.
## Successfully loaded ArchRProject!## Computing overlap...## Success!4.7 Generate regulons
A long format data frame, representing the inferred regulons, is then generated. The data frame consists of three columns:
- tf (transcription factor)
- target gene
- peak to gene correlation between tf and target gene
## DataFrame with 2776926 rows and 10 columns
## idxATAC chr start end idxRNA target corr distance idxTF tf
## <integer> <factor> <integer> <integer> <integer> <character> <matrix> <numeric> <integer> <character>
## 1 15 chr1 912762 913262 7 NOC2L 46297 4 AGO1
## 2 15 chr1 912762 913262 7 NOC2L 46297 11 ARID4B
## 3 15 chr1 912762 913262 7 NOC2L 46297 12 ARID5B
## 4 15 chr1 912762 913262 7 NOC2L 46297 30 BCOR
## 5 15 chr1 912762 913262 7 NOC2L 46297 36 BRD4
## ... ... ... ... ... ... ... ... ... ... ...
## 2776922 210665 chrX 154815200 154815700 23515 F8 211490 1146 NFRKB
## 2776923 210665 chrX 154815200 154815700 23515 F8 211490 1175 POLR2H
## 2776924 210665 chrX 154815200 154815700 23515 F8 211490 1273 ZBTB8A
## 2776925 210665 chrX 154815200 154815700 23515 F8 211490 1456 ZNF589
## 2776926 210665 chrX 154815200 154815700 23515 F8 211490 1457 ZNF5924.8 (Optional) Annotate with TF motifs
So far the gene regulatory network was constructed from TF ChIP-seq exclusively. Some users would prefer to further annotate the regulatory elements with the presence of motifs. If motif annotation has been previously performed by ArchR, addMotifScore can retrieve this annotation from the ArchR project.
If motifs are available for a factor and the RE contains a motif, we return 1. If motifs are available and the RE does not contain a motif, we return 0. If no motifs are known for this particular factor (as in the case of co-factors or chromatin modifiers), we return NA.
If the user has not performed motif annotation with ArchR, we can also annotate the peaks with motifs using the Cisbp database (default) or user-provided PWMS. See ?addMotifScore
It is important to note that filtering for the presence of motifs removes a large fraction of the target genes. Motifs are often present in a small subset of the ChIP-seq peaks (can be as low as 10%). Second, indirect TF binding, possibly through its interaction partners, may have a true biological function. In this example, we continue with regulons containing the motifs, regulon.motif. However, if the user prefers to retain all target genes including REs without the motifs, they should proceed with regulon.
## retrieving motif information from ArchR project4.9 Prune network
Epiregulon prunes the network by performing tests of independence on the observed number of cells jointly expressing transcription factor (TF), regulatory element (RE) and target gene (TG) vs the expected number of cells if TF/RE and TG are independently expressed. We implement two tests, the binomial test and the chi-square test. In the binomial test, the expected probability is P(TF, RE) * P(TG), and the number of trials is the total number of cells, and the observed successes is the number of cells jointly expressing all three elements. In the chi-square test, the expected probability for having all 3 elements active is also P(TF, RE) * P(TG) and the probability otherwise is 1- P(TF, RE) * P(TG). The observed cell count for the active category is the number of cells jointly expressing all three elements, and the cell count for the inactive category is n - n_triple.
We calculate cluster-specific p-values if users supply cluster labels. This is useful if we are interested in cluster-specific networks. The pruned regulons can then be used to visualize differential networks for transcription factors of interest. See section on differential networks.
4.10 Add Weights
While the `pruneRegulon’ function provides statistics on the joint occurrence of TF-RE-TG, we would like to further estimate the strength of regulation. Biologically, this can be interpreted as the magnitude of gene expression changes induced by transcription factor activity. Epiregulon estimates the regulatory potential using one of the three measures: 1) correlation between TG and TF or between TG and the product of TF and RE, 2) mutual information between TG and TF expression or between TG and the product of TF and RE, or 3) Wilcoxon test statistics of target gene expression in cells jointly expressing all 3 elements vs cells that do not.
Two measures (correlation and Wilcoxon) give both the magnitude and directionality of changes whereas mutational information is always positive. The correlation and mutual information statistics are computed on grouped pseudobulks by user-supplied cluster labels and yield a single weight across all clusters per each TF-RE-target triplet. In contrast, the Wilcoxon method group cells based on the joint expression of TF, RE and TG in each single cell or in cell aggregates. If cell labels are provided, we calculate cluster-specific weights in addition to estimating weights from all the cells. Cell aggregation uses a default value of 10 cells and can help overcome sparsity and speed up computation. If cluster labels are provided, we can obtain weights of individual clusters and all cells combined.
In the case of drug treatment, however, the activity of TF is suppressed often not by downregulation of the TF gene expression, but by direct interference of the TF protein function. In this dataset, the drug enzalutamide blocks the ligand binding domain of the androgen receptor and prevents it from binding to the chromatin. As a result, while the AR gene expression stays the same, the chromatin accessibility of AR, as computed by chromVar in the ArchR package, is greatly reduced by 48 hour treatment of enzalutamide.
First, we visualize the AR expression and show that enzalutamide does not decrease AR expression.
library(epiregulon.extra)
plotActivityDim(sce = GeneExpressionMatrix,
activity_matrix = assay(GeneExpressionMatrix),
tf = "AR",
dimtype = "UMAP_Combined",
label = "label",
point_size = 1,
legend.label = "gene expression")
Then we extract the chromVarMatrix from ArchR project and then visualize the chromatin accessibility at AR bound sites. We can see that 48 hour of enzalutamide treatment reduced chromatin accessibility at AR bound sites
chromVarMatrix <- getMatrixFromProject(
ArchRProj = proj,
useMatrix = "MotifMatrix",
useSeqnames = NULL,
verbose = TRUE,
binarize = FALSE,
threads = 1
)## 2024-03-19 08:20:45.565593 : Organizing colData, 0.512 mins elapsed.## 2024-03-19 08:20:45.73947 : Organizing rowData, 0.515 mins elapsed.## 2024-03-19 08:20:45.742636 : Organizing rowRanges, 0.515 mins elapsed.## 2024-03-19 08:20:45.748942 : Organizing Assays (1 of 2), 0.515 mins elapsed.## 2024-03-19 08:20:45.847062 : Organizing Assays (2 of 2), 0.517 mins elapsed.## 2024-03-19 08:20:45.943753 : Constructing SummarizedExperiment, 0.519 mins elapsed.## 2024-03-19 08:20:49.662956 : Finished Matrix Creation, 0.581 mins elapsed.plotActivityDim(sce = GeneExpressionMatrix,
activity_matrix = assay(chromVarMatrix, "z"),
tf = "AR_689",
dimtype = "UMAP_Combined",
label = "label",
point_size = 1,
legend.label = "chromVar")
Next, we are going to compare 3 different weight methods. In the first method, the wilcoxon test compares target gene expression in cells meeting both the TF expression and accessibility cutoffs vs cells failing either the TF expression or/and accessibility cutoffs. Next, we try out the correlation method which comes in two flavors. When tf_re.merge = FALSE, weight is computed on the correlation of target gene expression vs TF gene expression. When tf_re.merge = TRUE, weight is computed on the correlation of target gene expression vs the product of TF expression and chromatin accessibility at TF-bound regulatory elements.
regulon.w.wilcox <- addWeights(regulon = pruned.regulon,
expMatrix = GeneExpressionMatrix,
exp_assay = "counts",
peakMatrix = PeakMatrix,
peak_assay = "counts",
clusters = GeneExpressionMatrix$label,
method = "wilcoxon")## adding weights using wilcoxon...regulon.w.corr <- addWeights(regulon = pruned.regulon,
expMatrix = GeneExpressionMatrix,
exp_assay = "counts",
peakMatrix = PeakMatrix,
peak_assay = "counts",
clusters = GeneExpressionMatrix$label,
method = "corr")## adding weights using corr...## calculating average expression across clusters...## computing weights...regulon.w.corr.re <- addWeights(regulon = pruned.regulon,
expMatrix = GeneExpressionMatrix,
exp_assay = "counts",
peakMatrix = PeakMatrix,
peak_assay = "counts",
clusters = GeneExpressionMatrix$label,
method = "corr",
tf_re.merge = TRUE)## adding weights using corr...## calculating average expression across clusters...## computing weights...4.11 Calculate TF activity
Finally, the activities for a specific TF in each cell are computed by averaging the weighted expressions of target genes linked to the TF. \[y=\frac{1}{n}\sum_{i=1}^{n} x_i * weight_i\] where \(y\) is the activity of a TF for a cell \(n\) is the total number of targets for a TF \(x_i\) is the log count expression of target i where i in {1,2,…,n} \(weight_i\) is the weight of TF and target i
We calculate three different activities corresponding to the different weighted regulons
score.combine.wilcox <- calculateActivity(expMatrix = GeneExpressionMatrix,
exp_assay = "counts",
regulon = regulon.w.wilcox,
normalize = TRUE,
mode = "weight",
method = "weightedMean")## calculating TF activity from regulon using weightedmean## Warning in calculateActivity(expMatrix = GeneExpressionMatrix, exp_assay = "counts", : The weight column contains multiple
## subcolumns but no cluster information was provided. Using first column to compute activity...## aggregating regulons...## creating weight matrix...## calculating activity scores...## normalize by mean...## normalize by the number of targets...score.combine.corr <- calculateActivity(expMatrix = GeneExpressionMatrix,
exp_assay = "counts",
regulon = regulon.w.corr,
normalize = TRUE,
mode = "weight",
method = "weightedMean")## calculating TF activity from regulon using weightedmean## aggregating regulons...## creating weight matrix...## calculating activity scores...## normalize by mean...## normalize by the number of targets...score.combine.corr.re <- calculateActivity(expMatrix = GeneExpressionMatrix,
exp_assay = "counts",
regulon = regulon.w.corr.re,
normalize = TRUE,
mode = "weight",
method = "weightedMean")## calculating TF activity from regulon using weightedmean## aggregating regulons...## creating weight matrix...## calculating activity scores...## normalize by mean...## normalize by the number of targets...We visualize the different activities side by side.
library(epiregulon.extra)
plotActivityViolin(activity_matrix = score.combine.wilcox,
tf = c( "AR"),
clusters = GeneExpressionMatrix$label) + ggtitle ("AR activity by wilcoxon")
plotActivityViolin(activity_matrix = score.combine.corr,
tf = c( "AR"),
clusters = GeneExpressionMatrix$label) + ggtitle ("AR activity by corr TF vs TG")
plotActivityViolin(activity_matrix = score.combine.corr.re,
tf = c( "AR"),
clusters = GeneExpressionMatrix$label) + ggtitle ("AR activity by corr TF*RE vs TG")
In this case, activity calculated from correlation based on TF and TG expression is clearly wrong because we see increased AR activity after Enzalutamide treatment despite it being an AR antagonist. Therefore, for drug treatment which often decouples TF gene expression and its activity, it is important to take into consideration both TF gene expression and RE chromatin accessibility; the latter may be a better indicator of TF function if the TF has an effect on the chromatin accessibility. In this case, the recommended methods are either wilcox or corr with tf_re.merge = TRUE.
The astute users could however detect a difference in the prediction of the AR activity in the resistant clones “RES-A” and “RES-B” with respect to the parental “LNCaP” between the two methods. For example, the corr with tf_re.merge = TRUE shows increased AR activity in “RES-B” compared to “LNCaP” because “RES-B” shows increased AR expression. In contrast, the wilcoxon method did not predict an increase in AR activity in “RES-B” because “RES-B” still shows reduced chromatin accessibility compared to “LNCaP”. Since wilcoxon takes into account the co-occurrence of both TF gene expression and RE chromatin accessibility, this method does not predict an overall increase in AR activity.
In the absence of the ground truth, it is difficult to judge which method is superior. Therefore, it is always crucial to validate key findings with additional empirical evidence. The most important disclaimer we wish to make is that all predictions by epiregulon should be robustly tested experimentally.
4.12 Perform differential activity
For the remaining steps, we continue with activity derived from the wilcoxon method.
markers <- findDifferentialActivity(activity_matrix = score.combine.wilcox,
clusters = GeneExpressionMatrix$label,
pval.type = "some",
direction = "up",
test.type = "t")Take the top differential TFs
## Using a logFC cutoff of 0.1 for class LNCaP for direction equal to any## Using a logFC cutoff of 0 for class LNCaP–ENZ48 for direction equal to any## Using a logFC cutoff of 0 for class LNCaP RES-A for direction equal to any## Using a logFC cutoff of 0 for class LNCaP RES-B for direction equal to any4.13 Visualize the results
First visualize the known differential TFs by bubble plot
plotBubble(activity_matrix = score.combine.wilcox,
tf = c("AR","FOXA1", "MYC","JUN"),
clusters = GeneExpressionMatrix$label)
Then visualize the most differential TFs by clusters
plotBubble(activity_matrix = score.combine.wilcox,
tf = markers.sig$tf,
clusters = GeneExpressionMatrix$label)## Warning in max(logpval[is.finite(logpval)]): no non-missing arguments to max; returning -Inf## Warning: Removed 52 rows containing missing values or values outside the scale range (`geom_point()`).
Visualize the known differential TFs by UMAP
plotActivityDim(sce = GeneExpressionMatrix,
activity_matrix = score.combine.wilcox,
tf = c( "AR", "FOXA1", "MYC", "JUN"),
dimtype = "UMAP_Combined",
label = "label",
point_size = 1,
ncol = 2,
nrow = 2)
Visualize the newly discovered differential TFs by UMAP
plotActivityDim(sce = GeneExpressionMatrix,
activity_matrix = score.combine.wilcox,
tf = c("SPDEF","HES4","ATF5","NR2F2"),
dimtype = "UMAP_Combined",
label = "label",
point_size = 1,
ncol = 2,
nrow = 2)
Visualize regulons by heatmap
rowData(GeneExpressionMatrix) <- NULL
plotHeatmapRegulon(sce=GeneExpressionMatrix,
tfs= c( "AR", "FOXA1", "MYC", "JUN"),
regulon=regulon.w.wilcox,
regulon_cutoff=0.1,
downsample=1000,
cell_attributes="label",
col_gap="label",
exprs_values="counts",
name="regulon heatmap",
column_title_rot = 45)
plotHeatmapActivity(activity=score.combine.wilcox,
sce=GeneExpressionMatrix,
tfs=rownames(score.combine.wilcox),
downsample=1000,
cell_attributes="label",
col_gap="label",
name = "transcription factor activity",
column_title_rot = 45)
## Geneset enrichment
Sometimes we are interested to know what pathways are enriched in the regulon of a particular TF. We can perform geneset enrichment using the enricher function from clusterProfiler.
#retrieve genesets
H <- EnrichmentBrowser::getGenesets(org = "hsa",
db = "msigdb",
cat = "H",
gene.id.type = "SYMBOL",
cache = FALSE)
C6 <- EnrichmentBrowser::getGenesets(org = "hsa",
db = "msigdb",
cat = "C6",
gene.id.type = "SYMBOL",
cache = FALSE)
#combine genesets and convert genesets to be compatible with enricher
gs <- c(H,C6)
gs.list <- do.call(rbind,lapply(names(gs),
function(x) {data.frame(gs=x, genes=gs[[x]])}))
enrichresults <- regulonEnrich(TF = c("AR", "FOXA1", "MYC", "JUN"),
regulon = regulon.w.wilcox,
weight = "weight",
weight_cutoff = 0,
genesets = gs.list)## AR## FOXA1## MYC## JUN
We can visualize the genesets of known factors as a network
plotGseaNetwork(tf = names(enrichresults),
enrichresults = enrichresults,
p.adj_cutoff = 0.1,
ntop_pathways = 10)
We can visualize the genesets of differential factors as a network
enrichresults <- regulonEnrich(TF = markers.sig$tf,
regulon = regulon.w.wilcox,
weight = "weight",
weight_cutoff = 0,
genesets = gs.list)## HES4## NFIB## SPDEF## NFYB## CEBPG## ATF5## BPTF## SREBF1## REST## JUN
## JUN## NR2F2## SMARCC1## NR2F6## ATF5plotGseaNetwork(tf = names(enrichresults),
enrichresults = enrichresults,
p.adj_cutoff = 0.1,
ntop_pathways = 10)
4.14 Differential network analysis
In addition to looking at the summed TF activity, a second approach to investigate differential TF activity is to compare and contrast target genes or network topology. In this example, we know that AR is downregulated in the Enzalutamide treated cells compared to parental LNCaP.
plotDiffNetwork(regulon.w.wilcox,
cutoff = 0,
tf = c("AR"),
weight = "weight",
clusters = c("LNCaP","LNCaP–ENZ48"),
layout = "stress")## Building graph using weight as edge weights
We perform edge subtracted graph between two conditions and rank TFs by degree centrality. In this example, positive centrality indicates higher activity in parental LNCaP and negative centrality indicates higher activity in Enzalutamide treated cells.
# construct a graph of the parental and enzalutamide treated cells respectively
LNCaP_network <- buildGraph(regulon.w.wilcox, weights = "weight", cluster="LNCaP")## Building graph using weight as edge weights## Building graph using weight as edge weights# construct a difference graph
diff_graph <- buildDiffGraph(LNCaP_network, ENZ_network, abs_diff = FALSE)
diff_graph <- addCentrality(diff_graph)
diff_graph <- normalizeCentrality(diff_graph)
rank_table <- rankTfs(diff_graph)
library(ggplot2)
ggplot(rank_table, aes(x = rank, y = centrality)) +
geom_point() +
ggrepel::geom_text_repel(data = rbind(head(rank_table,5),
tail(rank_table,5)),
aes(label = tf),
nudge_x = 0, nudge_y = 0, box.padding = 0.5) +
theme_classic() + ggtitle ("differential TFs (LNCaP-ENZ) ranked by degree centrality")
Sometimes, we are interested to identify interaction partners of the TFs of interest. This can be achieved by comparing the overlap of the targets genes for all the TFs and identify the most similar TFs by Jaccard similarity. To illustrate this function, we take a look at the top most similar 20 TFs to AR.
library(igraph)
diff_graph_filter <- subgraph.edges(diff_graph,
E(diff_graph)[E(diff_graph)$weight>0],
del=TRUE)
# compute a similarity matrix of all TFs
similarity_score <- calculateJaccardSimilarity(diff_graph_filter)
# Focus on AR
similarity_score_AR <- similarity_score[, "AR"]
similarity_df <- data.frame(similarity = head(sort(similarity_score_AR,
decreasing = TRUE),20),
TF = names(head(sort(similarity_score_AR,
decreasing = TRUE),20)))
similarity_df$TF <- factor(similarity_df$TF,
levels = rev(unique(similarity_df$TF)))
# plot top TFs most similar to SPI1
topTFplot <- ggplot(similarity_df, aes(x=TF, y=similarity)) +
geom_bar(stat="identity") +
coord_flip() +
ggtitle("AR similarity") +
theme_classic()
print(topTFplot)
In order to convince ourselves that our differential network is statistically significant, we permute the edges and obtain a background graph from averaging many iterations. Here, we plot the differential network graph subtracted by permuted graphs.
# create a permuted graph by rewiring the edges 100 times
permute_matrix <- permuteGraph(diff_graph_filter, "AR", 100, p=1)
permute_matrix <- permute_matrix[names(similarity_score_AR),]
diff_matrix <- similarity_score_AR-rowMeans(permute_matrix)
diff_matrix_df <- data.frame(similarity = head(sort(diff_matrix,
decreasing = TRUE),20),
TF = names(head(sort(diff_matrix,
decreasing = TRUE),20)))
diff_matrix_df$TF <- factor(diff_matrix_df$TF, levels = rev(unique(diff_matrix_df$TF)))
# plot top TFs most similar to AR
topTFplot <- ggplot(diff_matrix_df, aes(x=TF, y=similarity)) +
geom_bar(stat="identity") +
coord_flip() +
ggtitle("background subtracted AR similarity ") +
theme_classic()
print(topTFplot)
# obtain empirical p-values
p_matrix <- rowMeans(apply(permute_matrix, 2, function(x) {x > similarity_score_AR}))
p_matrix[names(head(sort(diff_matrix,decreasing = TRUE),20))]## JUND MYC HOXB13 FOXA1 NFIC CEBPB XBP1 MAZ ATF4 CTCF REST GATA2 ETV1 FOXP1 CEBPG JUN YY1
## 0.00 0.00 0.00 0.01 0.03 0.01 0.00 0.03 0.01 0.00 0.01 0.01 0.03 0.01 0.01 0.03 0.00
## NFIB ZNF148 USF2
## 0.00 0.00 0.014.15 Session Info
## R Under development (unstable) (2023-12-04 r85659)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 18.04.6 LTS
##
## Matrix products: default
## BLAS: /usr/local/lib/R/lib/libRblas.so
## LAPACK: /usr/local/lib/R/lib/libRlapack.so; LAPACK version 3.11.0
##
## Random number generation:
## RNG: L'Ecuyer-CMRG
## Normal: Inversion
## Sample: Rejection
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Etc/UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] parallel stats4 grid stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] org.Hs.eg.db_3.18.0 AnnotationDbi_1.65.2 nabor_0.5.0
## [4] epiregulon.archr_0.99.2 igraph_2.0.3 msigdbr_7.5.1
## [7] epiregulon.extra_0.99.5 BSgenome.Hsapiens.UCSC.hg38_1.4.5 BSgenome_1.71.2
## [10] rtracklayer_1.63.1 BiocIO_1.13.0 Biostrings_2.71.4
## [13] XVector_0.43.1 epiregulon_0.99.5 scMultiome_1.3.0
## [16] SingleCellExperiment_1.25.0 MultiAssayExperiment_1.29.1 ExperimentHub_2.11.1
## [19] AnnotationHub_3.11.1 BiocFileCache_2.11.1 dbplyr_2.4.0
## [22] BiocStyle_2.31.0 rhdf5_2.47.6 SummarizedExperiment_1.33.3
## [25] Biobase_2.63.0 RcppArmadillo_0.12.8.1.0 Rcpp_1.0.12
## [28] Matrix_1.6-4 GenomicRanges_1.55.3 GenomeInfoDb_1.39.9
## [31] IRanges_2.37.1 S4Vectors_0.41.4 BiocGenerics_0.49.1
## [34] sparseMatrixStats_1.15.0 MatrixGenerics_1.15.0 matrixStats_1.2.0
## [37] data.table_1.15.2 stringr_1.5.1 plyr_1.8.9
## [40] magrittr_2.0.3 ggplot2_3.5.0 gtable_0.3.4
## [43] gtools_3.9.5 gridExtra_2.3 devtools_2.4.5
## [46] usethis_2.2.3 ArchR_1.0.3
##
## loaded via a namespace (and not attached):
## [1] fs_1.6.3 bitops_1.0-7 enrichplot_1.23.1 DirichletMultinomial_1.45.0
## [5] TFBSTools_1.41.0 doParallel_1.0.17 RColorBrewer_1.1-3 HDO.db_0.99.1
## [9] httr_1.4.7 Rgraphviz_2.47.0 profvis_0.3.8 tools_4.4.0
## [13] backports_1.4.1 utf8_1.2.4 R6_2.5.1 HDF5Array_1.31.6
## [17] lazyeval_0.2.2 GetoptLong_1.0.5 rhdf5filters_1.15.4 urlchecker_1.0.1
## [21] withr_3.0.0 cli_3.6.2 Cairo_1.6-2 scatterpie_0.2.1
## [25] labeling_0.4.3 sass_0.4.9 KEGGgraph_1.63.0 readr_2.1.5
## [29] yulab.utils_0.1.4 Rsamtools_2.19.3 gson_0.1.0 DOSE_3.29.2
## [33] R.utils_2.12.3 scater_1.31.2 sessioninfo_1.2.2 limma_3.59.6
## [37] rstudioapi_0.15.0 RSQLite_2.3.5 shape_1.4.6.1 gridGraphics_0.5-1
## [41] generics_0.1.3 dplyr_1.1.4 GO.db_3.18.0 ggbeeswarm_0.7.2
## [45] fansi_1.0.6 abind_1.4-5 R.methodsS3_1.8.2 lifecycle_1.0.4
## [49] yaml_2.3.8 edgeR_4.1.18 qvalue_2.35.0 SparseArray_1.3.4
## [53] blob_1.2.4 promises_1.2.1 dqrng_0.3.2 crayon_1.5.2
## [57] miniUI_0.1.1.1 lattice_0.22-5 beachmat_2.19.1 cowplot_1.1.3
## [61] annotate_1.81.2 KEGGREST_1.43.0 magick_2.8.3 ComplexHeatmap_2.19.0
## [65] pillar_1.9.0 knitr_1.45 metapod_1.11.1 fgsea_1.29.0
## [69] rjson_0.2.21 codetools_0.2-19 fastmatch_1.1-4 glue_1.7.0
## [73] ggfun_0.1.4 remotes_2.5.0 treeio_1.27.0 vctrs_0.6.5
## [77] png_0.1-8 poweRlaw_0.80.0 cachem_1.0.8 xfun_0.42
## [81] S4Arrays_1.3.6 mime_0.12 tidygraph_1.3.1 pracma_2.4.4
## [85] iterators_1.0.14 statmod_1.5.0 bluster_1.13.0 ellipsis_0.3.2
## [89] nlme_3.1-164 ggtree_3.11.1 bit64_4.0.5 filelock_1.0.3
## [93] bslib_0.6.1 irlba_2.3.5.1 vipor_0.4.7 colorspace_2.1-0
## [97] seqLogo_1.69.0 DBI_1.2.2 ggrastr_1.0.2 tidyselect_1.2.1
## [101] bit_4.0.5 compiler_4.4.0 curl_5.2.1 graph_1.81.0
## [105] BiocNeighbors_1.21.2 DelayedArray_0.29.9 shadowtext_0.1.3 bookdown_0.38
## [109] checkmate_2.3.1 scales_1.3.0 caTools_1.18.2 rappdirs_0.3.3
## [113] digest_0.6.35 motifmatchr_1.25.0 rmarkdown_2.26 htmltools_0.5.7
## [117] pkgconfig_2.0.3 highr_0.10 fastmap_1.1.1 GlobalOptions_0.1.2
## [121] rlang_1.1.3 htmlwidgets_1.6.4 shiny_1.8.0 DelayedMatrixStats_1.25.1
## [125] farver_2.1.1 jquerylib_0.1.4 jsonlite_1.8.8 BiocParallel_1.37.1
## [129] GOSemSim_2.29.1 R.oo_1.26.0 BiocSingular_1.19.0 RCurl_1.98-1.14
## [133] ggplotify_0.1.2 scuttle_1.13.1 GenomeInfoDbData_1.2.11 patchwork_1.2.0
## [137] Rhdf5lib_1.25.1 munsell_0.5.0 ape_5.7-1 babelgene_22.9
## [141] viridis_0.6.5 EnrichmentBrowser_2.33.1 stringi_1.8.3 ggraph_2.2.1
## [145] MASS_7.3-60.1 zlibbioc_1.49.3 pkgbuild_1.4.4 ggrepel_0.9.5
## [149] CNEr_1.39.0 graphlayouts_1.1.1 splines_4.4.0 circlize_0.4.16
## [153] hms_1.1.3 locfit_1.5-9.9 reshape2_1.4.4 ScaledMatrix_1.11.1
## [157] pkgload_1.3.4 TFMPvalue_0.0.9 BiocVersion_3.19.1 XML_3.99-0.16.1
## [161] evaluate_0.23 scran_1.31.3 BiocManager_1.30.22 foreach_1.5.2
## [165] tweenr_2.0.3 tzdb_0.4.0 httpuv_1.6.14 polyclip_1.10-6
## [169] tidyr_1.3.1 purrr_1.0.2 clue_0.3-65 ggforce_0.4.2
## [173] rsvd_1.0.5 xtable_1.8-4 restfulr_0.0.15 tidytree_0.4.6
## [177] later_1.3.2 viridisLite_0.4.2 tibble_3.2.1 aplot_0.2.2
## [181] clusterProfiler_4.11.0 memoise_2.0.1 beeswarm_0.4.0 GenomicAlignments_1.39.4
## [185] cluster_2.1.6 GSEABase_1.65.1