Help

At any time, remember you an access the help files for any function (i.e. ?cisTarget), and the other tutorials in the package with the following commands:

# Explore tutorials in the web browser:
browseVignettes(package="RcisTarget") 

# Commnad line-based:
vignette(package="RcisTarget") # list
vignette("RcisTarget") # open

Reports

To generate an HTML report with your own data and comments, you can use the .Rmd file of this tutorial as template (i.e. copy the .Rmd file, and edit it as R notebook in RStudio).

vignetteFile <- paste(file.path(system.file('doc', package='RcisTarget')), 
                      "RcisTarget.Rmd", sep="/")
# Copy to edit as markdown
file.copy(vignetteFile, ".")
# Alternative: extract R code
Stangle(vignetteFile)

Gene-motif rankings

The score of each pair of gene-motif can be performed using different parameters. Therefore, we provide multiple databases (motif-rankings, or alternative: mirror), according to several possibilities:

• Species: Species of the input gene set. Available values: Human (Homo sapiens), mouse (Mus musculus) or fly (Drosophila melanogaster)
• Scoring/search space: determine the search space around the transcription start site (TSS). Available values: 500 bp uptream the TSS, 5kbp or 10kbp around the TSS (e.g. 10kbp upstream and 10kbp downstream of the TSS).
• Motif collection: We recommend to use the latest one (mc10_clust).

If you don’t know which database to choose, for an analisis of a gene list we would suggest using the 500bp upstream TSS, and a larger search space (e.g. TSS+/-5kbp or TSS+/-10kbp). Of course, selecting Human, Mouse or fly depending on your input gene list.

For other settings (e.g. new species), you can check the tutorial on how to create databases.

Each database is stored in a .feather file. Note that the download size is typically over 1GB (100GB for human region databases), we recommend downloading the files with zsync_curl` (see the Help with downloads).

Once you have the .feather files, they can be loaded with importRankings():

# Search space: 10k bp around TSS - HUMAN
motifRankings <- importRankings("hg38_10kbp_up_10kbp_down_full_tx_v10_clust.genes_vs_motifs.rankings.feather")
# Load the annotation to human transcription factors
data(motifAnnotations_hgnc)

Motif annotations

All the calculations performed by RcisTarget are motif-based. However, most users will be interested on TFs potentially regulating the gene-set. The association of motif to transcription factors are provided in an independent file. In addition to the motifs annotated by the source datatabse (i.e. direct annotation), we have also inferred some further annotations based on motif similarity and gene ortology (e.g. similarity to other genes annotated to the motif). These annotations are typically used with the functions cisTarget() or addMotifAnnotation().

For the motifs in version ‘mc10nr’ of the rankings, these annotations are already included in RcisTarget package and can be loaded with the command:

# mouse:
# data(motifAnnotations_mgi)

# human:
data(motifAnnotations_hgnc)
motifAnnotations[199:202,]

## Key: <motif, TF>
##            motif     TF directAnnotation inferred_Orthology inferred_MotifSimil
##           <char> <char>           <lgcl>             <lgcl>              <lgcl>
## 1: cisbp__M00192  SMAD3            FALSE               TRUE               FALSE
## 2: cisbp__M00193   TCF7            FALSE               TRUE               FALSE
## 3: cisbp__M00194   HBP1            FALSE               TRUE               FALSE
## 4: cisbp__M00195    CIC            FALSE               TRUE               FALSE
##        annotationSource
##                  <fctr>
## 1: inferredBy_Orthology
## 2: inferredBy_Orthology
## 3: inferredBy_Orthology
## 4: inferredBy_Orthology
##                                                                                                         description
##                                                                                                              <char>
## 1: gene is orthologous to ENSMUSG00000032402 in M. musculus (identity = 100%) which is directly annotated for motif
## 2:  gene is orthologous to ENSMUSG00000000782 in M. musculus (identity = 93%) which is directly annotated for motif
## 3:  gene is orthologous to ENSMUSG00000002996 in M. musculus (identity = 95%) which is directly annotated for motif
## 4:                       motif is annotated for orthologous gene ENSMUSG00000005442 in M. musculus (identity = 92%)
##    keptInRanking
##           <lgcl>
## 1:         FALSE
## 2:         FALSE
## 3:         FALSE
## 4:          TRUE

For other versions of the rankings, the function importAnnotations allows to import the annotations from the source file.

These annotations can be easily queried to obtain further information about specific motifs or TFs:

showLogo(motifAnnotations[(directAnnotation==TRUE) & keptInRanking 
                               & (TF %in% c("HIF1A", "HIF2A", "EPAS1")), ])

# From Bioconductor
if (!requireNamespace("BiocManager", quietly=TRUE))
    install.packages("BiocManager")
BiocManager::install("RcisTarget.hg19.motifDBs.cisbpOnly.500bp")

library(RcisTarget.hg19.motifDBs.cisbpOnly.500bp)
# Rankings
data(hg19_500bpUpstream_motifRanking_cispbOnly)
motifRankings <- hg19_500bpUpstream_motifRanking_cispbOnly
motifRankings

## Rankings for RcisTarget.
##   Organism: human (hg19)
##   Number of genes: 22284 (22285 available in the full DB)
##   Number of MOTIF: 4687
## ** This database includes rankings up to 5050
## 
## Subset (4.6k cisbp motifs) of the Human database scoring motifs 500bp upstream the TSS (hg19-500bp-upstream-7species.mc9nr)

# Annotations
data(hg19_motifAnnotation_cisbpOnly) 
motifAnnotations <- hg19_motifAnnotation_cisbpOnly

1. Calculate enrichment

The first step to estimate the over-representation of each motif on the gene-set is to calculate the Area Under the Curve (AUC) for each pair of motif-geneSet. This is calculated based on the recovery curve of the gene-set on the motif ranking (genes ranked decreasingly by the score of motif in its proximity, as provided in the motifRanking database).

motifs_AUC <- calcAUC(geneLists, motifRankings, nCores=1)

The AUC is provided as a matrix of Motifs by GeneSets. In principle, the AUC is mostly meant as input for the next step. However, it is also possible to explore the distribution of the scores, for example in a gene-set of interest:

auc <- getAUC(motifs_AUC)["hypoxia",]
hist(auc, main="hypoxia", xlab="AUC histogram",
     breaks=100, col="#ff000050", border="darkred")
nes3 <- (3*sd(auc)) + mean(auc)
abline(v=nes3, col="red")

2. Select significant motifs and/or annotate to TFs

The selection of significant motifs is done based on the Normalized Enrichment Score (NES). The NES is calculated -for each motif- based on the AUC distribution of all the motifs for the gene-set [(x-mean)/sd]. Those motifs that pass the given threshold (3.0 by default) are considered significant.

Furthermore, this step also allows to add the TFs annotated to the motif.

motifEnrichmentTable <- addMotifAnnotation(motifs_AUC, nesThreshold=3,
     motifAnnot=motifAnnotations,
     highlightTFs=list(hypoxia=c("HIF1A", "EPAS1")))

class(motifEnrichmentTable)

## [1] "data.table" "data.frame"

dim(motifEnrichmentTable)

## [1] 17  8

head(motifEnrichmentTable[,-"TF_lowConf", with=FALSE])

##    geneSet        motif   NES    AUC highlightedTFs TFinDB
##     <char>       <char> <num>  <num>         <char> <char>
## 1: hypoxia cisbp__M6275  4.41 0.0966   HIF1A, EPAS1     **
## 2: hypoxia cisbp__M0062  3.57 0.0841   HIF1A, EPAS1       
## 3: hypoxia cisbp__M6279  3.56 0.0840   HIF1A, EPAS1       
## 4: hypoxia cisbp__M6212  3.49 0.0829   HIF1A, EPAS1     **
## 5: hypoxia cisbp__M2332  3.48 0.0828   HIF1A, EPAS1       
## 6: hypoxia cisbp__M0387  3.41 0.0817   HIF1A, EPAS1       
##                   TF_highConf
##                        <char>
## 1: HIF1A (directAnnotation). 
## 2:                           
## 3: HMGA1 (directAnnotation). 
## 4: EPAS1 (directAnnotation). 
## 5:                           
## 6:

The cathegories considered high/low confidence can me modified with the arguments motifAnnot_highConfCat and motifAnnot_lowConfCat.

3. Identify the genes with the best enrichment for each Motif

Since RcisTarget searches for enrichment of a motif within a gene list, finding a motif ‘enriched’ does not imply that all the genes in the gene-list have a high score for the motif. In this way, the third step of the workflow is to identify which genes (of the gene-set) are highly ranked for each of the significant motifs.

There are two methods to identify these genes: (1) equivalent to the ones used in iRegulon and i-cisTarget (method="iCisTarget", recommended if running time is not an issue), and (2) a faster implementation based on an approximate distribution using the average at each rank (method="aprox", useful to scan multiple gene sets).

IMPORTANT: Make sure that the motifRankings are the same as in Step 1.

motifEnrichmentTable_wGenes <- addSignificantGenes(motifEnrichmentTable,
                      rankings=motifRankings, 
                      geneSets=geneLists)
dim(motifEnrichmentTable_wGenes)

## [1] 17 11

Plot for a few motifs:

geneSetName <- "hypoxia"
selectedMotifs <- c(sample(motifEnrichmentTable$motif, 2))
par(mfrow=c(2,2))
getSignificantGenes(geneLists[[geneSetName]], 
                    rankings=motifRankings,
                    signifRankingNames=selectedMotifs,
                    plotCurve=TRUE, maxRank=5000, genesFormat="none",
                    method="aprox")