diff --git a/.github/workflows/check-bioc.yml b/.github/workflows/check-bioc.yml new file mode 100644 index 0000000..8eca8c2 --- /dev/null +++ b/.github/workflows/check-bioc.yml @@ -0,0 +1,13 @@ +name: R-CMD-check-bioc + +on: + push: + branches: + - main + - master + - devel + pull_request: + +jobs: + run-check: + uses: bioc-package-rescue/workflows/.github/workflows/check-bioc.yml@main diff --git a/DESCRIPTION b/DESCRIPTION index 62c149b..5196c6f 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -4,9 +4,12 @@ Title: RcisTarget Identify transcription factor binding motifs enriched on a list of genes or genomic regions Version: 1.23.1 Date: 2024-03-02 -Author: Sara Aibar, Gert Hulselmans, Stein Aerts. Laboratory of Computational Biology. - VIB-KU Leuven Center for Brain & Disease Research. Leuven, Belgium -Maintainer: Gert Hulselmans +Authors@R: c( + person("Sara", "Aibar", role = "aut"), + person("Gert", "Hulselmans", email = "Gert.Hulselmans@kuleuven.be", role = c("aut", "cre")), + person("Stein", "Aerts", role = "aut"), + person("Laboratory of Computational Biology. VIB-KU Leuven Center for Brain & Disease Research. Leuven, Belgium", role = "cph") + ) Description: RcisTarget identifies transcription factor binding motifs (TFBS) over-represented on a gene list. In a first step, RcisTarget selects DNA motifs that are significantly @@ -31,6 +34,7 @@ Imports: dplyr, tibble, GSEABase, + igraph, methods, R.utils, stats, @@ -48,7 +52,6 @@ Suggests: foreach, gplots, rtracklayer, - igraph, knitr, RcisTarget.hg19.motifDBs.cisbpOnly.500bp, rmarkdown, diff --git a/R/02_addMotifAnnotation.R b/R/02_addMotifAnnotation.R index b8e5fdb..246a6ce 100755 --- a/R/02_addMotifAnnotation.R +++ b/R/02_addMotifAnnotation.R @@ -186,7 +186,7 @@ addMotifAnnotation <- function(auc, nesThreshold=3.0, digits=3, if(!is.null(motifAnnot)) { motifAnnot_subset <- motifAnnot[(motifAnnot[[idColumn]] %in% aucTable[[idColumn]]) - & (motifAnnot$TF %in% TFs), ][,c(idColumn, "TF", "annotationSource"),with=F] + & (motifAnnot$TF %in% TFs), ][,c(idColumn, "TF", "annotationSource"),with=FALSE] motifAnnot_subset <- split(motifAnnot_subset, motifAnnot_subset[[idColumn]]) for(motifName in names(motifAnnot_subset)) { diff --git a/R/aux_addLogo.R b/R/aux_addLogo.R index 9b86327..6aca62d 100755 --- a/R/aux_addLogo.R +++ b/R/aux_addLogo.R @@ -94,7 +94,7 @@ showLogo <- function(motifEnrDT, stop("Please indicate the column containing the motif id (argument 'motifCol') or set it to NULL.") } # TODO - met <- met[grep("transfac_pro__", met[[motifCol]], invert = T),] + met <- met[grep("transfac_pro__", met[[motifCol]], invert = TRUE),] } # For numeric columns, show only the number of significant digits... @@ -106,7 +106,7 @@ showLogo <- function(motifEnrDT, if(!is.null(colsToShow)) { colsToShow <- unique(unname(unlist(colsToShow))) colsToShow <- colsToShow[which(colsToShow %in% colnames(met))] - met <- met[, colsToShow, with=F] + met <- met[, colsToShow, with=FALSE] } # Show... diff --git a/R/aux_convertToTargetRegions.R b/R/aux_convertToTargetRegions.R index bccd95a..9156252 100644 --- a/R/aux_convertToTargetRegions.R +++ b/R/aux_convertToTargetRegions.R @@ -14,11 +14,13 @@ #' \code{vignette("RcisTarget")} #' #' @examples -#' \dontrun{ -#' ## To apply on a list of regionSets: -#' regionSets_db <- lapply(regionSets, function(x) -#' convertToTargetRegions(queryRegions=x, targetRegions=dbRegionsLoc)) -#' } +#' # Create mock query and target regions +#' library(GenomicRanges) +#' queryRegions <- GRanges("chr1", IRanges(start = c(10, 100), end = c(50, 150))) +#' targetRegions <- GRanges("chr1", IRanges(start = c(5, 80, 200), end = c(45, 120, 250))) +#' +#' # Convert query regions to target regions +#' convertToTargetRegions(queryRegions = queryRegions, targetRegions = targetRegions) #' @rdname convertToTargetRegions #' @importFrom methods isClass #' @importFrom GenomeInfoDb keepSeqlevels diff --git a/R/aux_getDbRegionsLoc.R b/R/aux_getDbRegionsLoc.R index 4f1bf5b..d536437 100644 --- a/R/aux_getDbRegionsLoc.R +++ b/R/aux_getDbRegionsLoc.R @@ -15,10 +15,24 @@ #' \code{vignette("RcisTarget")} #' #' @examples -#' \dontrun{ -#' featherFilePath <- "~/databases/dm6-regions-11species.mc9nr.feather" -#' dbRegionsLoc <- getDbRegionsLoc(featherFilePath) -#' } +#' library(arrow) +#' # Create a dummy data frame with drosophila-style region column names +#' df <- data.frame( +#' features = c("motif1", "motif2"), +#' `motif1__chr3L:1000-2000` = c(1, 2), +#' `motif2__chr3L:3000-4000` = c(3, 4), +#' check.names = FALSE +#' ) +#' +#' # Write to a temporary feather file +#' tmpFile <- tempfile(fileext = ".feather") +#' write_feather(df, tmpFile) +#' +#' # Get the region locations +#' getDbRegionsLoc(tmpFile) +#' +#' # Clean up +#' file.remove(tmpFile) #' @rdname getDbRegionsLoc #' @importFrom GenomeInfoDb sortSeqlevels seqlevels #' @importFrom GenomicRanges GRanges diff --git a/R/priv_calcErn.R b/R/priv_calcErn.R index af5d455..1b70210 100755 --- a/R/priv_calcErn.R +++ b/R/priv_calcErn.R @@ -44,7 +44,7 @@ maxRankExtra <- maxRank+nMean gsRankings.asMat <- as.matrix(gsRankings) # Much faster! gsRankings.asMat[gsRankings.asMat>maxRankExtra] <- NA - globalMat <- matrix(0, nrow=max(nrow(gsRankings),max(gsRankings.asMat, na.rm=T)), ncol=maxRankExtra) # nrow=nrow(gsRankings): can be out of range, add extra rows and subset at the end + globalMat <- matrix(0, nrow=max(nrow(gsRankings),max(gsRankings.asMat, na.rm=TRUE)), ncol=maxRankExtra) # nrow=nrow(gsRankings): can be out of range, add extra rows and subset at the end # x <- x[seq_len(min(length(x), nrow(globalMat)))] # or: # if(nrow(globalMat) < length(x)) x <- x[seq_len(nrow(globalMat))] # cannot be in the loop... too slow! for(i in 1:nrow(gsRankings)) # (TO DO: Paralellize?) diff --git a/inst/examples/example_addSignificantGenes.R b/inst/examples/example_addSignificantGenes.R index 0ed27f5..932ecd0 100755 --- a/inst/examples/example_addSignificantGenes.R +++ b/inst/examples/example_addSignificantGenes.R @@ -61,16 +61,17 @@ motifEnr_wIncidMat <- addSignificantGenes(motifEnrichmentTable, genesFormat = "incidMatrix") motifEnr_wIncidMat <- as.data.frame(motifEnr_wIncidMat) -which(colnames(motifEnr_wIncidMat) == "rankAtMax") +rankAtMaxIdx <- match("rankAtMax", colnames(motifEnr_wIncidMat)) +if(is.na(rankAtMaxIdx)) stop("Column 'rankAtMax' not found in results.") -incidMat <- motifEnr_wIncidMat[,8:ncol(motifEnr_wIncidMat)] +incidMat <- as.matrix(motifEnr_wIncidMat[, (rankAtMaxIdx + 1):ncol(motifEnr_wIncidMat)]) rownames(incidMat) <- motifEnr_wIncidMat[,"motif"] incidMat <- incidMat[, colSums(incidMat)>0, drop=FALSE] # Plot as network par(mfrow=c(1,1)) library(igraph) -plot(graph.incidence(incidMat)) +plot(graph_from_incidence_matrix(incidMat)) ############################################################### # Alternative method: getSignificantGenes() diff --git a/inst/examples/example_workflow.R b/inst/examples/example_workflow.R index cce10fb..0d559c2 100755 --- a/inst/examples/example_workflow.R +++ b/inst/examples/example_workflow.R @@ -1,38 +1,33 @@ # RcisTarget workflow for advanced users: # Running the workflow steps individually -\dontrun{ - -################################################## -#### Load your gene sets # As example, the package includes an Hypoxia gene set: txtFile <- paste(file.path(system.file('examples', package='RcisTarget')), "hypoxiaGeneSet.txt", sep="/") geneLists <- list(hypoxia=read.table(txtFile, stringsAsFactors=FALSE)[,1]) - -#### Load databases -## Motif rankings: Select according to organism and distance around TSS -## (See the vignette for URLs to download) -motifRankings <- importRankings("~/databases/hg38_10kbp_up_10kbp_down_full_tx_v10_clust.genes_vs_motifs.rankings.feather") -## Motif - TF annotation: -data("motifAnnotations_hgnc") # human TFs (for motif collection 10) -################################################## - -#### Run RcisTarget +# Since the suggested database package is in Suggests, +# we guard the example to avoid hard failures if it is not installed: +if (requireNamespace("RcisTarget.hg19.motifDBs.cisbpOnly.500bp", quietly = TRUE)) { + + # Load rankings database + data(hg19_500bpUpstream_motifRanking_cispbOnly, package="RcisTarget.hg19.motifDBs.cisbpOnly.500bp") + motifRankings <- hg19_500bpUpstream_motifRanking_cispbOnly -# Step 1. Calculate AUC -motifs_AUC <- calcAUC(geneLists, motifRankings) + # Load motif annotations + data(motifAnnotations_hgnc_v9) # human TFs (for motif collection 9) + motifAnnotation <- motifAnnotations_hgnc_v9 -# Step 2. Select significant motifs, add TF annotation & format as table -motifEnrichmentTable <- addMotifAnnotation(motifs_AUC, - motifAnnot=motifAnnotations) + # Step 1. Calculate AUC + motifs_AUC <- calcAUC(geneLists, motifRankings) -# Step 3 (optional). Identify genes that have the motif significantly enriched -# (i.e. genes from the gene set in the top of the ranking) -motifEnrichmentTable_wGenes <- addSignificantGenes(motifEnrichmentTable, - geneSets=geneLists, - rankings=motifRankings, - method="aprox") + # Step 2. Select significant motifs, add TF annotation & format as table + motifEnrichmentTable <- addMotifAnnotation(motifs_AUC, + motifAnnot=motifAnnotation) + # Step 3 (optional). Identify genes that have the motif significantly enriched + motifEnrichmentTable_wGenes <- addSignificantGenes(motifEnrichmentTable, + geneSets=geneLists, + rankings=motifRankings, + method="aprox") } diff --git a/man/addSignificantGenes.Rd b/man/addSignificantGenes.Rd index efebbbe..55e7c97 100755 --- a/man/addSignificantGenes.Rd +++ b/man/addSignificantGenes.Rd @@ -255,9 +255,10 @@ geneLists <- list(hypoxia=read.table(txtFile, stringsAsFactors=FALSE)[,1]) ## (See the vignette for URLs to download) # motifRankings <- importRankings("hg19-500bp-upstream-7species.mc9nr.feather") +if (requireNamespace("RcisTarget.hg19.motifDBs.cisbpOnly.500bp", quietly = TRUE)) { + ## For this example we will use a SUBSET of the ranking/motif databases: -library(RcisTarget.hg19.motifDBs.cisbpOnly.500bp) -data(hg19_500bpUpstream_motifRanking_cispbOnly) +data(hg19_500bpUpstream_motifRanking_cispbOnly, package="RcisTarget.hg19.motifDBs.cisbpOnly.500bp") motifRankings <- hg19_500bpUpstream_motifRanking_cispbOnly ## Motif - TF annotation: @@ -302,16 +303,17 @@ motifEnr_wIncidMat <- addSignificantGenes(motifEnrichmentTable, genesFormat = "incidMatrix") motifEnr_wIncidMat <- as.data.frame(motifEnr_wIncidMat) -which(colnames(motifEnr_wIncidMat) == "rankAtMax") +rankAtMaxIdx <- match("rankAtMax", colnames(motifEnr_wIncidMat)) +if(is.na(rankAtMaxIdx)) stop("Column 'rankAtMax' not found in results.") -incidMat <- motifEnr_wIncidMat[,8:ncol(motifEnr_wIncidMat)] +incidMat <- as.matrix(motifEnr_wIncidMat[, (rankAtMaxIdx + 1):ncol(motifEnr_wIncidMat)]) rownames(incidMat) <- motifEnr_wIncidMat[,"motif"] incidMat <- incidMat[, colSums(incidMat)>0, drop=FALSE] # Plot as network par(mfrow=c(1,1)) library(igraph) -plot(graph.incidence(incidMat)) +plot(graph_from_incidence_matrix(incidMat)) ############################################################### # Alternative method: getSignificantGenes() @@ -322,8 +324,7 @@ onlyGenes <- getSignificantGenes(geneSet=geneLists$hypoxia, plotCurve=TRUE, rankings=motifRankings, method="aprox") - - +} } \seealso{ Previous step in the workflow: \code{\link{addMotifAnnotation}}. diff --git a/man/calcAUC.Rd b/man/calcAUC.Rd index 3c96694..fb2be10 100755 --- a/man/calcAUC.Rd +++ b/man/calcAUC.Rd @@ -107,41 +107,36 @@ that are significantly over-represented in the gene-set. # RcisTarget workflow for advanced users: # Running the workflow steps individually -\dontrun{ - -################################################## -#### Load your gene sets # As example, the package includes an Hypoxia gene set: txtFile <- paste(file.path(system.file('examples', package='RcisTarget')), "hypoxiaGeneSet.txt", sep="/") geneLists <- list(hypoxia=read.table(txtFile, stringsAsFactors=FALSE)[,1]) - -#### Load databases -## Motif rankings: Select according to organism and distance around TSS -## (See the vignette for URLs to download) -motifRankings <- importRankings("hg19-500bp-upstream-7species.mc9nr.feather") - -## Motif - TF annotation: -data(motifAnnotations_hgnc_v9) # human TFs (for motif collection 9) -motifAnnotation <- motifAnnotations_hgnc_v9 -################################################## - -#### Run RcisTarget - -# Step 1. Calculate AUC -motifs_AUC <- calcAUC(geneLists, motifRankings) - -# Step 2. Select significant motifs, add TF annotation & format as table -motifEnrichmentTable <- addMotifAnnotation(motifs_AUC, - motifAnnot=motifAnnotation) - -# Step 3 (optional). Identify genes that have the motif significantly enriched -# (i.e. genes from the gene set in the top of the ranking) -motifEnrichmentTable_wGenes <- addSignificantGenes(motifEnrichmentTable, - geneSets=geneLists, - rankings=motifRankings, - method="aprox") +# Since the suggested database package is in Suggests, +# we guard the example to avoid hard failures if it is not installed: +if (requireNamespace("RcisTarget.hg19.motifDBs.cisbpOnly.500bp", quietly = TRUE)) { + + # Load rankings database + data(hg19_500bpUpstream_motifRanking_cispbOnly, package="RcisTarget.hg19.motifDBs.cisbpOnly.500bp") + motifRankings <- hg19_500bpUpstream_motifRanking_cispbOnly + + # Load motif annotations + data(motifAnnotations_hgnc_v9) # human TFs (for motif collection 9) + motifAnnotation <- motifAnnotations_hgnc_v9 + + # Step 1. Calculate AUC + motifs_AUC <- calcAUC(geneLists, motifRankings) + + # Step 2. Select significant motifs, add TF annotation & format as table + motifEnrichmentTable <- addMotifAnnotation(motifs_AUC, + motifAnnot=motifAnnotation) + + # Step 3 (optional). Identify genes that have the motif significantly enriched + # (i.e. genes from the gene set in the top of the ranking) + motifEnrichmentTable_wGenes <- addSignificantGenes(motifEnrichmentTable, + geneSets=geneLists, + rankings=motifRankings, + method="aprox") } } \seealso{ diff --git a/man/convertToTargetRegions.Rd b/man/convertToTargetRegions.Rd index 91e5a7d..342edbb 100644 --- a/man/convertToTargetRegions.Rd +++ b/man/convertToTargetRegions.Rd @@ -34,11 +34,13 @@ IDs of the regions in the "target regions" overlapping with the query regions. Convert a set of input regions to the overlapping regions in the target set. } \examples{ -\dontrun{ - ## To apply on a list of regionSets: - regionSets_db <- lapply(regionSets, function(x) - convertToTargetRegions(queryRegions=x, targetRegions=dbRegionsLoc)) - } +# Create mock query and target regions +library(GenomicRanges) +queryRegions <- GRanges("chr1", IRanges(start = c(10, 100), end = c(50, 150))) +targetRegions <- GRanges("chr1", IRanges(start = c(5, 80, 200), end = c(45, 120, 250))) + +# Convert query regions to target regions +convertToTargetRegions(queryRegions = queryRegions, targetRegions = targetRegions) } \seealso{ \code{\link{getDbRegionsLoc}}. diff --git a/man/getDbRegionsLoc.Rd b/man/getDbRegionsLoc.Rd index a740966..45c932a 100644 --- a/man/getDbRegionsLoc.Rd +++ b/man/getDbRegionsLoc.Rd @@ -28,8 +28,22 @@ For \bold{human/mouse} the region locations are stored in a separate object (i.e. \code{}). This function is not needed. } \examples{ -\dontrun{ -featherFilePath <- "~/databases/dm6-regions-11species.mc9nr.feather" -dbRegionsLoc <- getDbRegionsLoc(featherFilePath) -} +library(arrow) +# Create a dummy data frame with drosophila-style region column names +df <- data.frame( + features = c("motif1", "motif2"), + `motif1__chr3L:1000-2000` = c(1, 2), + `motif2__chr3L:3000-4000` = c(3, 4), + check.names = FALSE +) + +# Write to a temporary feather file +tmpFile <- tempfile(fileext = ".feather") +write_feather(df, tmpFile) + +# Get the region locations +getDbRegionsLoc(tmpFile) + +# Clean up +file.remove(tmpFile) }