Contents

1 Introduction

Gene regulatory networks model the underlying gene regulation hierarchies that drive gene expression and cell states. The main functions of this package are to construct gene regulatory networks and infer transcription factor (TF) activity at the single cell level by integrating scATAC-seq and scRNA-seq data and incorporating of public bulk TF ChIP-seq data.

There are three related packages: epiregulon, epiregulon.extra and epiregulon.archr, the two of which are available through Bioconductor and the last of which is only available through github. The basic epiregulon package takes in gene expression and chromatin accessibility matrices in the form of SingleCellExperiment objects, constructs gene regulatory networks (“regulons”) and outputs the activity of transcription factors at the single cell level. The epiregulon.extra package provides a suite of tools for enrichment analysis of target genes, visualization of target genes and transcription factor activity, and network analysis which can be run on the epiregulon output. If the users would like to start from ArchR projects instead of SingleCellExperiment objects, they may choose to use epiregulon.archr package, which allows for seamless integration with the ArchR package, and continue with the rest of the workflow offered in epiregulon.extra.

For full documentation of the epiregulon package, please refer to the epiregulon book.

2 Installation

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

Alternatively, you could install from github

devtools::install_github(repo ='xiaosaiyao/epiregulon')

Load package

library(epiregulon)
#> Loading required package: SingleCellExperiment
#> Loading required package: SummarizedExperiment
#> Loading required package: MatrixGenerics
#> Loading required package: matrixStats
#> 
#> Attaching package: 'MatrixGenerics'
#> The following objects are masked from 'package:matrixStats':
#> 
#>     colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
#>     colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
#>     colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
#>     colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
#>     colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
#>     colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
#>     colWeightedMeans, colWeightedMedians, colWeightedSds,
#>     colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
#>     rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
#>     rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
#>     rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
#>     rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
#>     rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
#>     rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
#>     rowWeightedSds, rowWeightedVars
#> Loading required package: GenomicRanges
#> Loading required package: stats4
#> Loading required package: BiocGenerics
#> 
#> Attaching package: 'BiocGenerics'
#> The following objects are masked from 'package:stats':
#> 
#>     IQR, mad, sd, var, xtabs
#> The following objects are masked from 'package:base':
#> 
#>     Filter, Find, Map, Position, Reduce, anyDuplicated, aperm, append,
#>     as.data.frame, basename, cbind, colnames, dirname, do.call,
#>     duplicated, eval, evalq, get, grep, grepl, intersect, is.unsorted,
#>     lapply, mapply, match, mget, order, paste, pmax, pmax.int, pmin,
#>     pmin.int, rank, rbind, rownames, sapply, setdiff, table, tapply,
#>     union, unique, unsplit, which.max, which.min
#> Loading required package: S4Vectors
#> 
#> Attaching package: 'S4Vectors'
#> The following object is masked from 'package:utils':
#> 
#>     findMatches
#> The following objects are masked from 'package:base':
#> 
#>     I, expand.grid, unname
#> Loading required package: IRanges
#> Loading required package: GenomeInfoDb
#> Loading required package: Biobase
#> Welcome to Bioconductor
#> 
#>     Vignettes contain introductory material; view with
#>     'browseVignettes()'. To cite Bioconductor, see
#>     'citation("Biobase")', and for packages 'citation("pkgname")'.
#> 
#> Attaching package: 'Biobase'
#> The following object is masked from 'package:MatrixGenerics':
#> 
#>     rowMedians
#> The following objects are masked from 'package:matrixStats':
#> 
#>     anyMissing, rowMedians

3 Data preparation

Prior to using epiregulon, single cell preprocessing needs to performed by user’s favorite methods. The following components are required:
1. Peak matrix from scATAC-seq containing the chromatin accessibility information
2. Gene expression matrix from either paired or unpaired scRNA-seq. RNA-seq integration needs to be performed for unpaired dataset.
3. Dimensionality reduction matrix from either single modality dataset or joint scRNA-seq and scATAC-seq

This tutorial demonstrates the basic functions of epiregulon, using the reprogram-seq dataset which can be downloaded from the scMultiome package. In this example, prostate cancer cells (LNCaP) were infected in separate wells with viruses encoding 4 transcription factors (NKX2-1, GATA6, FOXA1 and FOXA2) and a positive control (mNeonGreen) before pooling. The identity of the infected transcription factors was tracked through cell hashing (available in the field hash_assignment of the colData) and serves as the ground truth.

# load the MAE object
library(scMultiome)
#> Loading required package: AnnotationHub
#> Loading required package: BiocFileCache
#> Loading required package: dbplyr
#> 
#> Attaching package: 'AnnotationHub'
#> The following object is masked from 'package:Biobase':
#> 
#>     cache
#> Loading required package: ExperimentHub
#> Loading required package: MultiAssayExperiment

mae <- scMultiome::reprogramSeq()
#> see ?scMultiome and browseVignettes('scMultiome') for documentation
#> loading from cache

# extract peak matrix
PeakMatrix <- mae[["PeakMatrix"]]

# extract expression matrix
GeneExpressionMatrix <- mae[["GeneExpressionMatrix"]]
rownames(GeneExpressionMatrix) <- rowData(GeneExpressionMatrix)$name

# define the order of hash_assigment
GeneExpressionMatrix$hash_assignment <- 
  factor(as.character(GeneExpressionMatrix$hash_assignment),
         levels = c("HTO10_GATA6_UTR", "HTO2_GATA6_v2", "HTO8_NKX2.1_UTR", "HTO3_NKX2.1_v2",
                    "HTO1_FOXA2_v2", "HTO4_mFOXA1_v2", "HTO6_hFOXA1_UTR", "HTO5_NeonG_v2"))
# extract dimensional reduction matrix
reducedDimMatrix <- reducedDim(mae[['TileMatrix500']], "LSI_ATAC")

# transfer UMAP_combined from TileMatrix to GeneExpressionMatrix
reducedDim(GeneExpressionMatrix, "UMAP_Combined") <- 
  reducedDim(mae[['TileMatrix500']], "UMAP_Combined")

Visualize singleCellExperiment by UMAP


scater::plotReducedDim(GeneExpressionMatrix, 
                       dimred = "UMAP_Combined", 
                       text_by = "Clusters", 
                       colour_by = "Clusters")

4 Quick start

4.1 Retrieve bulk TF ChIP-seq binding sites

First, we retrieve a GRangesList object containing the binding sites of all the transcription factors and co-regulators. These binding sites are derived from bulk ChIP-seq data in the ChIP-Atlas and ENCODE databases. For the same transcription factor, multiple ChIP-seq files from different cell lines or tissues are merged. For further information on how these peaks are derived, please refer to ?epiregulon::getTFMotifInfo. Currently, human genomes hg19 and hg38 and mouse mm10 are supported.

grl <- getTFMotifInfo(genome = "hg38")
#> see ?scMultiome and browseVignettes('scMultiome') for documentation
#> loading from cache
grl
#> 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>

We recommend the use of ChIP-seq data over motif for estimating TF occupancy. However, if the user would like to start from motifs, it is possible to switch to the motif mode. The user can provide the peak regions as a GRanges object and the location of the motifs will be annotated based on Cisbp from chromVARmotifs (human_pwms_v2 and mouse_pwms_v2, version 0.2)

grl.motif <- getTFMotifInfo(genome = "hg38",  
                            mode = "motif", 
                            peaks = rowRanges(PeakMatrix))

4.3 Add TF motif binding to peaks

The next step is to add the TF binding information by overlapping regions of the peak matrix with the bulk chip-seq database. The output is a data frame object with three columns:

  1. idxATAC - index of the peak in the peak matrix
  2. idxTF - index in the gene expression matrix corresponding to the transcription factor
  3. tf - name of the transcription factor

overlap <- addTFMotifInfo(grl = grl, p2g = p2g, 
                          peakMatrix = PeakMatrix)
#> Computing overlap...
#> Success!
head(overlap)
#>      idxATAC idxTF     tf
#> 9312      37     4   AGO1
#> 9313      37     7 ARID1B
#> 9314      37     8  ARID2
#> 9315      37    18   ATF3
#> 9316      37    21   ATF7
#> 9317      37    23  BACH1

4.4 Generate regulons

A DataFrame, representing the inferred regulons, is then generated. The DataFrame consists of ten columns:

  1. idxATAC - index of the peak in the peak matrix
  2. chr - chromosome number
  3. start - start position of the peak
  4. end - end position of the peak
  5. idxRNA - index in the gene expression matrix corresponding to the target gene
  6. target - name of the target gene
  7. distance - distance between the transcription start site of the target gene and the middle of the peak
  8. idxTF - index in the gene expression matrix corresponding to the transcription factor
  9. tf - name of the transcription factor
  10. corr - correlation between target gene expression and the chromatin accessibility at the peak. if cluster labels are provided, this field is a matrix with columns names corresponding to correlation across all cells and for each of the clusters.

regulon <- getRegulon(p2g = p2g, overlap = overlap, 
                      aggregate = FALSE)
regulon
#> DataFrame with 743656 rows and 10 columns
#>          idxATAC         chr     start       end    idxRNA      target
#>        <integer> <character> <integer> <integer> <integer> <character>
#> 1             37        chr1   1020509   1021009        22  AL645608.4
#> 2             37        chr1   1020509   1021009        22  AL645608.4
#> 3             37        chr1   1020509   1021009        22  AL645608.4
#> 4             37        chr1   1020509   1021009        22  AL645608.4
#> 5             37        chr1   1020509   1021009        22  AL645608.4
#> ...          ...         ...       ...       ...       ...         ...
#> 743652    126541        chrX 154374945 154375445     36389        FLNA
#> 743653    126541        chrX 154374945 154375445     36389        FLNA
#> 743654    126590        chrX 155228844 155229344     36426       CLIC2
#> 743655    126590        chrX 155228844 155229344     36426       CLIC2
#> 743656    126590        chrX 155228844 155229344     36426       CLIC2
#>         distance     idxTF          tf     corr
#>        <integer> <integer> <character> <matrix>
#> 1         106338         4        AGO1 0.592745
#> 2         106338         7      ARID1B 0.592745
#> 3         106338         8       ARID2 0.592745
#> 4         106338        18        ATF3 0.592745
#> 5         106338        21        ATF7 0.592745
#> ...          ...       ...         ...      ...
#> 743652      3723      1492      ZNF692 0.572236
#> 743653      3723      1557        ZXDB 0.572236
#> 743654    105268       114       FOXA1 0.527470
#> 743655    105268       128       GATA2 0.527470
#> 743656    105268       952       SUMO2 0.527470

4.6 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 TF and target gene expression, 2) mutual information between the TF and target gene expression and 3) Wilcoxon test statistics of target gene expression in cells jointly expressing all 3 elements vs cells that do not.

Two of the measures (correlation and Wilcoxon statistics) give both the magnitude and directionality of changes whereas mutual information is always positive. The correlation and mutual information statistics are computed on pseudobulks aggregated by user-supplied cluster labels, whereas the Wilcoxon method groups cells into two categories: 1) the active category of cells jointly expressing TF, RE and TG and 2) the inactive category consisting of the remaining cells.

We calculate cluster-specific weights if users supply cluster labels.


regulon.w <- addWeights(regulon = pruned.regulon,
                        expMatrix  = GeneExpressionMatrix,
                        exp_assay  = "normalizedCounts",
                        peakMatrix = PeakMatrix,
                        peak_assay = "counts",
                        clusters = GeneExpressionMatrix$Clusters,
                        method = "wilcox")

regulon.w

4.7 (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 regulatory elements with the presence of motifs. We provide an option to annotate peaks with motifs from the Cisbp database. If no motifs are present for this particular factor (as in the case of co-factors or chromatin modifiers), we return NA. 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. The users can also provide their own motif annotation through the pwms argument.


regulon.w.motif <- addMotifScore(regulon = regulon.w,
                                 peaks = rowRanges(PeakMatrix),
                                 species = "human",
                                 genome = "hg38")
#> annotating peaks with motifs
#> see ?scMultiome and browseVignettes('scMultiome') for documentation
#> loading from cache
#> 
#> 
#> Attaching package: 'Biostrings'
#> The following object is masked from 'package:base':
#> 
#>     strsplit
#> 
#> Attaching package: 'rtracklayer'
#> The following object is masked from 'package:BiocIO':
#> 
#>     FileForFormat
#> The following object is masked from 'package:AnnotationHub':
#> 
#>     hubUrl

# if desired, set weight to 0 if no motif is found
regulon.w.motif$weight[regulon.w.motif$motif == 0] <- 0

regulon.w.motif
#> DataFrame with 3290 rows and 15 columns
#>        idxATAC         chr     start       end    idxRNA      target  distance
#>      <integer> <character> <integer> <integer> <integer> <character> <integer>
#> 1           86        chr1   1290083   1290583        58        DVL1     58527
#> 2          288        chr1   2526365   2526865       108  AL139246.1     33107
#> 3          640        chr1   7988538   7989038       208     TNFRSF9     47672
#> 4          654        chr1   8040536   8041036       212      ERRFI1     24903
#> 5          787        chr1   9290024   9290524       229       SPSB1      2368
#> ...        ...         ...       ...       ...       ...         ...       ...
#> 3286    120296        chr9  36118633  36119133     34311        CLTA     71721
#> 3287    121835        chr9  98192900  98193400     34702      CORO2A     20016
#> 3288    122254        chr9 108991440 108991940     34784     CTNNAL1     27885
#> 3289    122461        chr9 113174265 113174765     34825     SLC31A1     46777
#> 3290    123380        chr9 128635245 128635745     35034      SPTAN1     82687
#>          idxTF          tf     corr                                   pval
#>      <integer> <character> <matrix>                               <matrix>
#> 1          490          AR 0.524540 1.91265e-03:0.50585953:0.000249278:...
#> 2          490          AR 0.602193 5.02384e-02:0.00249757:1.000000000:...
#> 3          490          AR 0.875387 7.53206e-09:0.15249790:1.000000000:...
#> 4          490          AR 0.550906 9.18676e-02:0.76535176:1.000000000:...
#> 5          490          AR 0.891738 9.64012e-13:0.51358560:1.000000000:...
#> ...        ...         ...      ...                                    ...
#> 3286       807      NKX2-1 0.522644                    0.035831331:1:1:...
#> 3287       807      NKX2-1 0.541195                    0.000419607:1:1:...
#> 3288       807      NKX2-1 0.872888                    0.013099814:1:1:...
#> 3289       807      NKX2-1 0.683659                    0.045120208:1:1:...
#> 3290       807      NKX2-1 0.745530                    0.045120208:1:1:...
#>                               stats                qval            weight
#>                            <matrix>            <matrix>          <matrix>
#> 1     9.63153:0.4426222:13.4176:... 1.00000e+00:1:1:...         0:0:0:...
#> 2     3.83348:9.1423696: 0.0000:... 1.00000e+00:1:1:...         0:0:0:...
#> 3    33.39232:2.0470852: 0.0000:... 1.24761e-04:1:1:...         0:0:0:...
#> 4     2.84134:0.0890786: 0.0000:... 1.00000e+00:1:1:...         0:0:0:...
#> 5    50.91607:0.4267549: 0.0000:... 1.60836e-08:1:1:...         0:0:0:...
#> ...                             ...                 ...               ...
#> 3286                4.40511:0:0:...           1:1:1:... 0.0343753:0:0:...
#> 3287               12.44280:0:0:...           1:1:1:... 0.0000000:0:0:...
#> 3288                6.15558:0:0:...           1:1:1:... 0.0000000:0:0:...
#> 3289                4.01414:0:0:...           1:1:1:... 0.0000000:0:0:...
#> 3290                4.01414:0:0:...           1:1:1:... 0.0000000:0:0:...
#>          motif
#>      <numeric>
#> 1            0
#> 2            0
#> 3            0
#> 4            0
#> 5            0
#> ...        ...
#> 3286         1
#> 3287         0
#> 3288         0
#> 3289         0
#> 3290         0

4.8 (Optional) Annotate with log fold changes

It is sometimes helpful to filter the regulons based on gene expression changes between two conditions. The addLogFC function is a wrapper around scran::findMarkers and adds extra columns of log changes to the regulon DataFrame. The users can specify the reference condition in logFC_ref and conditions for comparison in logFC_condition. If these are not provided, log fold changes are calculated for every condition in the cluster labels against the rest of the conditions.

regulon.w <- addLogFC(regulon = regulon.w,
                      clusters = GeneExpressionMatrix$hash_assignment,
                      expMatrix  = GeneExpressionMatrix,
                      assay.type  = "normalizedCounts",
                      pval.type = "any",
                      logFC_condition = c("HTO10_GATA6_UTR", 
                                          "HTO2_GATA6_v2", 
                                          "HTO8_NKX2.1_UTR"),
                      logFC_ref = "HTO5_NeonG_v2")

regulon.w
#> DataFrame with 3290 rows and 23 columns
#>        idxATAC         chr     start       end    idxRNA      target  distance
#>      <integer> <character> <integer> <integer> <integer> <character> <integer>
#> 1           86        chr1   1290083   1290583        58        DVL1     58527
#> 2          288        chr1   2526365   2526865       108  AL139246.1     33107
#> 3          640        chr1   7988538   7989038       208     TNFRSF9     47672
#> 4          654        chr1   8040536   8041036       212      ERRFI1     24903
#> 5          787        chr1   9290024   9290524       229       SPSB1      2368
#> ...        ...         ...       ...       ...       ...         ...       ...
#> 3286    120296        chr9  36118633  36119133     34311        CLTA     71721
#> 3287    121835        chr9  98192900  98193400     34702      CORO2A     20016
#> 3288    122254        chr9 108991440 108991940     34784     CTNNAL1     27885
#> 3289    122461        chr9 113174265 113174765     34825     SLC31A1     46777
#> 3290    123380        chr9 128635245 128635745     35034      SPTAN1     82687
#>          idxTF          tf     corr                                   pval
#>      <integer> <character> <matrix>                               <matrix>
#> 1          490          AR 0.524540 1.91265e-03:0.50585953:0.000249278:...
#> 2          490          AR 0.602193 5.02384e-02:0.00249757:1.000000000:...
#> 3          490          AR 0.875387 7.53206e-09:0.15249790:1.000000000:...
#> 4          490          AR 0.550906 9.18676e-02:0.76535176:1.000000000:...
#> 5          490          AR 0.891738 9.64012e-13:0.51358560:1.000000000:...
#> ...        ...         ...      ...                                    ...
#> 3286       807      NKX2-1 0.522644                    0.035831331:1:1:...
#> 3287       807      NKX2-1 0.541195                    0.000419607:1:1:...
#> 3288       807      NKX2-1 0.872888                    0.013099814:1:1:...
#> 3289       807      NKX2-1 0.683659                    0.045120208:1:1:...
#> 3290       807      NKX2-1 0.745530                    0.045120208:1:1:...
#>                               stats                qval
#>                            <matrix>            <matrix>
#> 1     9.63153:0.4426222:13.4176:... 1.00000e+00:1:1:...
#> 2     3.83348:9.1423696: 0.0000:... 1.00000e+00:1:1:...
#> 3    33.39232:2.0470852: 0.0000:... 1.24761e-04:1:1:...
#> 4     2.84134:0.0890786: 0.0000:... 1.00000e+00:1:1:...
#> 5    50.91607:0.4267549: 0.0000:... 1.60836e-08:1:1:...
#> ...                             ...                 ...
#> 3286                4.40511:0:0:...           1:1:1:...
#> 3287               12.44280:0:0:...           1:1:1:...
#> 3288                6.15558:0:0:...           1:1:1:...
#> 3289                4.01414:0:0:...           1:1:1:...
#> 3290                4.01414:0:0:...           1:1:1:...
#>                                weight HTO10_GATA6_UTR.vs.HTO5_NeonG_v2.FDR
#>                              <matrix>                            <numeric>
#> 1    0.0509796:0.0300186:0.531885:...                          1.00000e+00
#> 2    0.0328174:0.1507898:0.000000:...                          8.87712e-01
#> 3    0.0903869:0.0673971:0.000000:...                          9.15316e-07
#> 4    0.0289190:0.0166923:0.000000:...                          8.83663e-06
#> 5    0.1169358:0.0323406:0.000000:...                          9.02310e-18
#> ...                               ...                                  ...
#> 3286                0.0343753:0:0:...                          1.00000e+00
#> 3287                0.0546077:0:0:...                          1.05779e-01
#> 3288                0.0471468:0:0:...                          2.66833e-19
#> 3289                0.0333575:0:0:...                          7.48668e-02
#> 3290                0.0277381:0:0:...                          8.87712e-01
#>      HTO10_GATA6_UTR.vs.HTO5_NeonG_v2.p.value
#>                                     <numeric>
#> 1                                   -0.428753
#> 2                                   -1.145259
#> 3                                  -10.189167
#> 4                                   -9.057294
#> 5                                  -22.259745
#> ...                                       ...
#> 3286                                -0.257374
#> 3287                                -3.833972
#> 3288                               -23.941330
#> 3289                                -4.063803
#> 3290                                -1.898997
#>      HTO10_GATA6_UTR.vs.HTO5_NeonG_v2.logFC HTO2_GATA6_v2.vs.HTO5_NeonG_v2.FDR
#>                                   <numeric>                          <numeric>
#> 1                                -0.0259827                        1.00000e+00
#> 2                                 0.0047951                        1.00000e+00
#> 3                                 0.2606361                        7.53057e-04
#> 4                                 0.7002420                        5.79312e-09
#> 5                                 0.8283787                        3.15502e-17
#> ...                                     ...                                ...
#> 3286                               0.036667                        8.48287e-01
#> 3287                              -0.156505                        3.24179e-02
#> 3288                               1.023267                        8.72255e-12
#> 3289                              -0.312983                        5.15569e-01
#> 3290                               0.205053                        3.07454e-03
#>      HTO2_GATA6_v2.vs.HTO5_NeonG_v2.p.value
#>                                   <numeric>
#> 1                                -0.0799202
#> 2                                 0.0000000
#> 3                                -6.8308339
#> 4                               -12.7390401
#> 5                               -21.7211619
#> ...                                     ...
#> 3286                               -2.15348
#> 3287                               -4.65304
#> 3288                              -15.86899
#> 3289                               -2.69505
#> 3290                               -6.02625
#>      HTO2_GATA6_v2.vs.HTO5_NeonG_v2.logFC HTO8_NKX2.1_UTR.vs.HTO5_NeonG_v2.FDR
#>                                 <numeric>                            <numeric>
#> 1                              0.00598935                             1.000000
#> 2                              0.00000000                             1.000000
#> 3                              0.11468090                             0.958744
#> 4                              0.94159934                             0.399423
#> 5                              0.71132340                             1.000000
#> ...                                   ...                                  ...
#> 3286                            -0.171173                          0.958743982
#> 3287                            -0.164758                          1.000000000
#> 3288                             0.975554                          0.000295046
#> 3289                            -0.211080                          0.958743982
#> 3290                             0.451261                          1.000000000
#>      HTO8_NKX2.1_UTR.vs.HTO5_NeonG_v2.p.value
#>                                     <numeric>
#> 1                                   -0.396194
#> 2                                    0.000000
#> 3                                   -1.146313
#> 4                                   -3.458271
#> 5                                   -0.108267
#> ...                                       ...
#> 3286                                -1.715298
#> 3287                                -0.808155
#> 3288                                -8.158592
#> 3289                                -1.647813
#> 3290                                -1.069575
#>      HTO8_NKX2.1_UTR.vs.HTO5_NeonG_v2.logFC
#>                                   <numeric>
#> 1                                -0.0238713
#> 2                                 0.0000000
#> 3                                 0.0068971
#> 4                                 0.2421805
#> 5                                -0.0053878
#> ...                                     ...
#> 3286                              0.1548462
#> 3287                             -0.0545147
#> 3288                              0.3115049
#> 3289                              0.1947646
#> 3290                              0.1178727

4.9 Calculate TF activity

Finally, the activities for a specific TF in each cell are computed by averaging expressions of target genes linked to the TF weighted by the test statistics of choice, chosen from either correlation, mutual information or the Wilcoxon test statistics. \[y=\frac{1}{n}\sum_{i=1}^{n} x_i * weights_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} and \(weights_i\) is the weight of TF - target \(i\)

score.combine <- calculateActivity(expMatrix = GeneExpressionMatrix,
                                   regulon = regulon.w, 
                                   mode = "weight", 
                                   method = "weightedMean", 
                                   exp_assay = "normalizedCounts",
                                   normalize = FALSE)
#> calculating TF activity from regulon using weightedmean
#> Warning in calculateActivity(expMatrix = GeneExpressionMatrix, regulon =
#> regulon.w, : 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 the number of targets...

score.combine[1:5,1:5]
#> 5 x 5 sparse Matrix of class "dgCMatrix"
#>        reprogram#TTAGGAACAAGGTACG-1 reprogram#GAGCGGTCAACCTGGT-1
#> AR                      0.043408741                   0.06010198
#> FOXA1                   0.035260181                   0.04408534
#> FOXA2                   0.008868738                   0.04479888
#> GATA6                   0.033885727                   0.13969050
#> NKX2-1                  0.065298024                   0.02791934
#>        reprogram#TTATAGCCACCCTCAC-1 reprogram#TGGTGATTCCTGTTCA-1
#> AR                       0.02734512                   0.03422265
#> FOXA1                    0.01708987                   0.01827346
#> FOXA2                    0.01153884                   0.01998275
#> GATA6                    0.01667011                   0.01904236
#> NKX2-1                   0.01410221                   0.01023911
#>        reprogram#TCGGTTCTCACTAGGT-1
#> AR                       0.04019227
#> FOXA1                    0.03426071
#> FOXA2                    0.01368200
#> GATA6                    0.02905760
#> NKX2-1                   0.05365239

5 Session Info

sessionInfo()
#> R version 4.4.0 (2024-04-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 22.04.4 LTS
#> 
#> Matrix products: default
#> BLAS:   /home/biocbuild/bbs-3.19-bioc/R/lib/libRblas.so 
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
#> 
#> 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       
#> 
#> time zone: America/New_York
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats4    stats     graphics  grDevices utils     datasets  methods  
#> [8] base     
#> 
#> other attached packages:
#>  [1] BSgenome.Hsapiens.UCSC.hg38_1.4.5 BSgenome_1.72.0                  
#>  [3] rtracklayer_1.64.0                BiocIO_1.14.0                    
#>  [5] Biostrings_2.72.1                 XVector_0.44.0                   
#>  [7] scMultiome_1.4.0                  MultiAssayExperiment_1.30.2      
#>  [9] ExperimentHub_2.12.0              AnnotationHub_3.12.0             
#> [11] BiocFileCache_2.12.0              dbplyr_2.5.0                     
#> [13] epiregulon_1.0.1                  SingleCellExperiment_1.26.0      
#> [15] SummarizedExperiment_1.34.0       Biobase_2.64.0                   
#> [17] GenomicRanges_1.56.1              GenomeInfoDb_1.40.1              
#> [19] IRanges_2.38.0                    S4Vectors_0.42.0                 
#> [21] BiocGenerics_0.50.0               MatrixGenerics_1.16.0            
#> [23] matrixStats_1.3.0                 BiocStyle_2.32.0                 
#> 
#> loaded via a namespace (and not attached):
#>   [1] bitops_1.0-7                filelock_1.0.3             
#>   [3] R.oo_1.26.0                 tibble_3.2.1               
#>   [5] XML_3.99-0.16.1             DirichletMultinomial_1.46.0
#>   [7] lifecycle_1.0.4             pwalign_1.0.0              
#>   [9] edgeR_4.2.0                 lattice_0.22-6             
#>  [11] backports_1.5.0             magrittr_2.0.3             
#>  [13] limma_3.60.2                sass_0.4.9                 
#>  [15] rmarkdown_2.27              jquerylib_0.1.4            
#>  [17] yaml_2.3.8                  metapod_1.12.0             
#>  [19] cowplot_1.1.3               DBI_1.2.3                  
#>  [21] CNEr_1.40.0                 abind_1.4-5                
#>  [23] zlibbioc_1.50.0             R.utils_2.12.3             
#>  [25] purrr_1.0.2                 RCurl_1.98-1.14            
#>  [27] pracma_2.4.4                rappdirs_0.3.3             
#>  [29] GenomeInfoDbData_1.2.12     ggrepel_0.9.5              
#>  [31] irlba_2.3.5.1               seqLogo_1.70.0             
#>  [33] annotate_1.82.0             dqrng_0.4.1                
#>  [35] DelayedMatrixStats_1.26.0   codetools_0.2-20           
#>  [37] DelayedArray_0.30.1         scuttle_1.14.0             
#>  [39] tidyselect_1.2.1            UCSC.utils_1.0.0           
#>  [41] farver_2.1.2                ScaledMatrix_1.12.0        
#>  [43] viridis_0.6.5               GenomicAlignments_1.40.0   
#>  [45] jsonlite_1.8.8              BiocNeighbors_1.22.0       
#>  [47] motifmatchr_1.26.0          scater_1.32.0              
#>  [49] tools_4.4.0                 TFMPvalue_0.0.9            
#>  [51] Rcpp_1.0.12                 glue_1.7.0                 
#>  [53] gridExtra_2.3               SparseArray_1.4.8          
#>  [55] xfun_0.44                   dplyr_1.1.4                
#>  [57] HDF5Array_1.32.0            withr_3.0.0                
#>  [59] BiocManager_1.30.23         fastmap_1.2.0              
#>  [61] rhdf5filters_1.16.0         bluster_1.14.0             
#>  [63] fansi_1.0.6                 caTools_1.18.2             
#>  [65] digest_0.6.35               rsvd_1.0.5                 
#>  [67] R6_2.5.1                    mime_0.12                  
#>  [69] colorspace_2.1-0            GO.db_3.19.1               
#>  [71] gtools_3.9.5                poweRlaw_0.80.0            
#>  [73] RSQLite_2.3.7               R.methodsS3_1.8.2          
#>  [75] utf8_1.2.4                  generics_0.1.3             
#>  [77] httr_1.4.7                  S4Arrays_1.4.1             
#>  [79] TFBSTools_1.42.0            pkgconfig_2.0.3            
#>  [81] gtable_0.3.5                blob_1.2.4                 
#>  [83] htmltools_0.5.8.1           bookdown_0.39              
#>  [85] scales_1.3.0                png_0.1-8                  
#>  [87] scran_1.32.0                knitr_1.47                 
#>  [89] tzdb_0.4.0                  reshape2_1.4.4             
#>  [91] rjson_0.2.21                checkmate_2.3.1            
#>  [93] curl_5.2.1                  cachem_1.1.0               
#>  [95] rhdf5_2.48.0                stringr_1.5.1              
#>  [97] BiocVersion_3.19.1          parallel_4.4.0             
#>  [99] vipor_0.4.7                 AnnotationDbi_1.66.0       
#> [101] restfulr_0.0.15             pillar_1.9.0               
#> [103] grid_4.4.0                  vctrs_0.6.5                
#> [105] BiocSingular_1.20.0         beachmat_2.20.0            
#> [107] xtable_1.8-4                cluster_2.1.6              
#> [109] beeswarm_0.4.0              evaluate_0.24.0            
#> [111] readr_2.1.5                 tinytex_0.51               
#> [113] magick_2.8.3                cli_3.6.2                  
#> [115] locfit_1.5-9.9              compiler_4.4.0             
#> [117] Rsamtools_2.20.0            rlang_1.1.4                
#> [119] crayon_1.5.2                labeling_0.4.3             
#> [121] plyr_1.8.9                  ggbeeswarm_0.7.2           
#> [123] stringi_1.8.4               viridisLite_0.4.2          
#> [125] BiocParallel_1.38.0         munsell_0.5.1              
#> [127] Matrix_1.7-0                hms_1.1.3                  
#> [129] sparseMatrixStats_1.16.0    bit64_4.0.5                
#> [131] ggplot2_3.5.1               Rhdf5lib_1.26.0            
#> [133] KEGGREST_1.44.0             statmod_1.5.0              
#> [135] highr_0.11                  igraph_2.0.3               
#> [137] memoise_2.0.1               bslib_0.7.0                
#> [139] bit_4.0.5