R version: R version 4.2.1 (2022-06-23)
Bioconductor version: 3.16
Package version: 1.12.0
Biological information encoded in the DNA sequence is often organized into independent modules or clusters. For instance, the eukaryotic system of expression relies on combinations of homotypic or heterotypic transcription factors (TFs) which play an important role in activating and repressing target genes. Identifying clusters of genomic features is a frequent task and includes application such as genomic regions enriched for the presence of a combination of transcription factor binding sites or enriched for the presence of mutations often referred as regions with increase mutational hotspots.
fcScan is designed to detect clusters of genomic features, based on user defined search criteria. Such criteria include:
fcScan is designed to handle large number of genomic features (i.e data generated by High Throughput Sequencing).
fcScan depends on the following packages:
Currently, fcScan has one main function, the getCluster function. Additional
functionality will be added for future releases including cross-species
identification of orthologous clusters.
The input for getCluster is given through the parameter x.
This function accepts input in data frame, GRanges as well as vector of files
in BED and VCF (compressed or not) formats. BED and VCF files are loaded
by packages rtracklayer and VariantAnnotation respectively. There is no
limit to the number of files the user can define.
When input is data frame or GRanges objects are given, they should contain the following “named” 5 columns:
| seqnames | start | end | strand | site |
|---|---|---|---|---|
| chr1 | 10 | 15 | + | a |
| chr1 | 17 | 20 | + | b |
| chr1 | 25 | 30 | + | b |
| chr1 | 27 | 35 | + | a |
| chr1 | 32 | 40 | + | c |
| chr1 | 41 | 48 | + | b |
The seqnames, start and end columns define the name of the chromosome, the starting and ending position of the feature in the chromosome, respectively. The strand column defines the strand. The site column contains the ID of the site and will be used for clustering. The start and end columns are numeric while the remaining columns are characters.
Note: when input is data frame, the data is considered zero-based
as in BED format.
Window size is set using w. It defines the maximum size of
the clusters.
The clustering condition c defines the number and name of genomic features to
search for based on features names defined in the site column
c = c("a" = 1, "b" = 2)
This searches for clusters with 3 sites, One a site and two b sites.
Another way of writing the condition, only if input is a vector to file paths,
is the following x = ("a.bed", "b.bed"), c = c(1,2)
Given 2 files, a.bed and b.bed, this condition states that the user is
looking for clusters made from 1 “a” site and 2 “b” sites. In this case, the
order of sites defined in c is relative to the order of files.
When input is a data frame or GRanges object (instead of files), the
user needs to explicitly define the site names along with the desired
number relative to each site.
For instance, giving the condition as c = c(1,2) for a data frame or GRanges
is not allowed.
x = dataframe, c = c("a" = 1, "b" = 2) where a and b are valid site
names in site column in the dataframe/GRanges
Users can exclude clusters containing a specific site(s).
This is done by specifying zero 0 in the condition as
c = c("a" = 1, "b" = 2, "c" = 0). In this case, any
cluster(s) containing c site will be excluded even if it has 1 a and 2 b
sites.
By default, clustering will be performed on both strands and on all seqnames
unless specified by the user using the s and seqnames arguments to limit
the search on a specific strand and/or seqname.
Users can choose to cluster on one specific seqname (seqnames = "chr1"), or
on several seqnames (seqnames = c("chr1","chr3","chr4"))
(Default for seqnames is NULL) meaning that clustering on all seqnames
will be performed.
For s, the values allowed are:
The gap/overlap between adjacent clusters, and not sites, can be controled
using the overlap option.
When overlap is a positive integer, adjacent clusters will be separated by
a minimum of the given value.
When overlap is negative, adjacent clusters will be allowed to overlap by
a maximum of the given value. (Default is set to 0)
greedy allows the user to control the number of genomic features found
in clusters.
When greedy = FALSE, getCluster will build clusters having the required
window size and will label TRUE the ones that contain the exact number
of sites provided in the condition argument.
Clusters having the user defined window size but not satisfying
the condition will be labelled as FALSE.
When greedy = TRUE, additional sites defined in condition will be added
to the cluster as long as the cluster size is below the defined window size.
(Default is set to FALSE)
The order option defines the desired order of the sites within identified
clusters. For instance, order = c("a","b","b")
will search for clusters containing 1 a and 2 b sites and checks if
they are in the specified order. Clusters with 1 a and 2 b sites that
do not contain the specified order will be rejected.
When greedy is set to TRUE, order can be satisfied if a subcluster contains
the desired order. For example if a cluster has a, a, b, b, b sites, it
satisfies the required order (a, b, b) and therefore will be considered as
a correct cluster . (Default is set to NULL)
The site_orientation option defines the orientation or strandness of sites
in the found clusters. This option cannot be used if order is NULL.
site_orientation should be specified for each site in order.
For instance, if order = c("a","b","b"), we can define site_orientation
for each site respectively as follow: site_orientation = c("+","-","-").
The cluster will be correct if it satisfies the required order and sites
orientation. (Default is set to NULL)
The site_overlap option defines the maximum or minimum distance allowed
between sites in the cluster.
When site_overlap is a positive integer, sites within a cluster should have
minimum distance and above, then the cluster is considered TRUE cluster.
When site_overlap is a negative integer, sites within a cluster should
have max distance and below, then the cluster is considered TRUE cluster.
When site_overlap is zero, distance between sites is not taken into
consideration. (Default is set to 0)
The cluster_by option allows the usage of different ways in scanning for
sites. Using this option, scanning does not ‘jump’ clusters found.
cluster_by can be startsEnds, endsStarts, starts, ends,
or middles. If we choose startsEnds, the clustering begins at the start of
the first site and it ends at the end of the last site within the window size.
If we choose endsStarts, scanning begins at the end of the
first site and it ends at the start of the last site within the window size.
if we choose starts, the clustering begins at the start of the first site
and it ends at the start of the last site within the window size.
if we choose ends, the clustering begins at the end of the first site
and it ends at the end of the last site within the window size.
Also, if we choose middles, the clustering begins at the middle
of the first site and it ends at the middle of the last site within the
window size. (Default is set to startsEnds)
The allow_clusters_overlap allows the user whether to include overlapping
clusters or not. (Default is set to FALSE)
This option can be only used in case where the “cluster_by” choosen is
“startsEnds” or “ends”, otherwise it is ignored. If it is set to TRUE,
partially overlapping sites within window size (at the end of the cluster)
are allowed. However, if it is set to FALSE, partially overlapping sites
are not included in the cluster found. (Default is set to FALSE)
This option allows the user to choose which partially overlapping site are
included. The amount of overlap between the last site found and the window
size controls the addition or removal of those sites. For example, if
partial_overlap_percentage is set to 0.6, this means that the window size
should overlap this site by 60% of its size. However, if it fails, this site
won’t be included in the cluster found. This option can be only used in case
where the “cluster_by” chosen is “startsEnds” or “ends”,
otherwise it is ignored. (Default is set to NULL)
This option allows to distribute the run over multiple threads. It is
recommended to use multiple threads when the input dataset is large.
When threads is set to 0, getCluster reserves the adequate number of
threads and performs the run. This will dramatically increase the speed of
getCluster. The user can choose the number of threads between 1 and the
maximum number of available threads. (Default is set to 1)
The verbose option allows the printing of additional messages and the
list of clusters that failed for lack of correct combination of sites.
This option is used mainly for debugging purposes. (Default is set to FALSE)
The output of getCluster is a GRanges object with fields:
The algorithm returns all clusters containing the correct count of
sites/features, unless verbose is set to TRUE. If the combination,
overlap and order options are satisfied, the cluster is considered a
TRUE cluster.
The status of a cluster can be either PASS, excludedSites, orderFail,
siteOrientation or siteOverlap
PASS is a cluster that satisfied the desired combination, overlap, order
and sites orientation.
orderFail is a cluster that had the required combination but did not
satisfy the required order of sites.
excludedSites is a cluster that had the required combination and order
but it has one or more sites to exclude.
siteOrientation is a cluster that had the required combination and order
but it has one or more sites with different orientation than requested.
siteOverlap is a cluster that meets all the requirements
but the distance between its sites is not respected.
NOTE: If the user is using greedy = FALSE and order contains values more
than in the condition parameter (c), an error will be raised.
However, if greedy = TRUE, then using order with more values than the
condition parameter is allowed since the cluster may contain more
sites than the required c condition as long as the window size is satisfied.
Example using getCluster:
getCluster looks for desired genomic regions of interest like
transcription factor binding sites, within a window size
and specific condition.
This function accepts a data frame and GRanges object. getCluster also accepts
BED or VCF (or mix of both) files as input.
The output of getCluster is a GRanges object that
contains the genomic coordinates (seqnames, ranges and strand) and
three metadata columns:
sites: contains clusters of sites that conforms with the condition c
specified in getCluster.
isCluster: TRUE if the cluster found conform with the condition c
and the order (if indicated in condition) and FALSE
if the cluster fails to conform with the condition or order.
status: PASS if isCluster equals TRUE. However, if isCluster is
FALSE, status shows why the found cluster is not a TRUE cluster. If the
order of sites is not respected in the found cluster, status would return
OrderFail. in Addition, if the cluster found contains non desired sites,
it returns excludedSites. Moreover, if the sites orientation is not
respected in found cluster, status would return siteOrientation.
Finally , if the distance between sites in the cluster found does not conform
to the condition specified by the user, status would return siteOverlap.
In this example, we ask getCluster to look for clusters that contains
one site “s1”, one site “s2” and zero “s3” sites. In addition, we requested
clusters to have sites in the order s1,s2 and having orientation
“+”,“+” respectively with minimum distance between sites of 2bp.
x1 = data.frame(seqnames = rep("chr1", times = 17),
start = c(1,10,17,25,27,32,41,47,60,70,87,94,99,107,113,121,132),
end = c(8,15,20,30,35,40,48,55,68,75,93,100,105,113,120,130,135),
strand = c("+","+","+","+","+","+","+","+","+",
"+","+","+","+","+","+","+","-"),
site = c("s3","s1","s2","s2","s1","s2","s1","s1","s2","s1","s2",
"s2","s1","s2","s1","s1","s2"))
clusters = getCluster(x1, w = 20, c = c("s1" = 1, "s2" = 1, "s3" = 0),
greedy = TRUE, order = c("s1","s2"), site_orientation=c("+","+"),
site_overlap = 2, verbose = TRUE)
#> 17 entries loaded
#> Running getCluster using 1 threads
#> Time difference of 0.1900337 secs
clusters
#> GRanges object with 9 ranges and 3 metadata columns:
#> seqnames ranges strand | sites isCluster status
#> <Rle> <IRanges> <Rle> | <character> <logical> <character>
#> [1] chr1 2-20 * | s3,s1,s2 FALSE excludedSites
#> [2] chr1 11-30 * | s1,s2,s2 TRUE PASS
#> [3] chr1 33-48 * | s2,s1 FALSE orderFail
#> [4] chr1 48-68 * | s1,s2 TRUE PASS
#> [5] chr1 88-105 * | s2,s2,s1 FALSE orderFail
#> [6] chr1 95-113 * | s2,s1,s2 FALSE siteOverlap
#> [7] chr1 100-120 * | s1,s2,s1 FALSE siteOverlap
#> [8] chr1 108-120 * | s2,s1 FALSE orderFail
#> [9] chr1 122-135 * | s1,s2 FALSE siteOrientation
#> -------
#> seqinfo: 1 sequence from an unspecified genome; no seqlengthsAnother example but using GRanges as input:
in this example, we ask getCluster to look for clusters that contains
one site s1 and two sites s2 within a window size of 25 bp. Also,
we requested clusters to be searched as + strand.
suppressMessages(library(GenomicRanges))
x = GRanges(
seqnames = Rle("chr1", 16),
ranges = IRanges(start = c(10L,17L,25L,27L,32L,41L,47L,
60L,70L,87L,94L,99L,107L,113L,121L,132L),
end = c(15L,20L,30L,35L,40L,48L,55L,68L,75L,93L,100L,105L,
113L,120L,130L,135L)),
strand = Rle("+",16),
site = c("s1","s2","s2","s1","s2","s1","s1","s2",
"s1","s2","s2","s1","s2","s1","s1","s2"))
clusters = getCluster(x, w = 25, c = c("s1"=1,"s2"=2), s = "+")
#> 16 entries loaded
#> Running getCluster using 1 threads
#> Time difference of 0.1090674 secs
clusters
#> GRanges object with 2 ranges and 3 metadata columns:
#> seqnames ranges strand | sites isCluster status
#> <Rle> <IRanges> <Rle> | <character> <logical> <character>
#> [1] chr1 10-30 + | s1,s2,s2 TRUE PASS
#> [2] chr1 87-105 + | s2,s2,s1 TRUE PASS
#> -------
#> seqinfo: 1 sequence from an unspecified genome; no seqlengthssessionInfo()
#> R version 4.2.1 (2022-06-23)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Ubuntu 20.04.5 LTS
#>
#> Matrix products: default
#> BLAS: /home/biocbuild/bbs-3.16-bioc/R/lib/libRblas.so
#> LAPACK: /home/biocbuild/bbs-3.16-bioc/R/lib/libRlapack.so
#>
#> 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
#>
#> attached base packages:
#> [1] stats4 stats graphics grDevices utils datasets methods
#> [8] base
#>
#> other attached packages:
#> [1] GenomicRanges_1.50.0 GenomeInfoDb_1.34.0 IRanges_2.32.0
#> [4] S4Vectors_0.36.0 BiocGenerics_0.44.0 fcScan_1.12.0
#> [7] BiocStyle_2.26.0
#>
#> loaded via a namespace (and not attached):
#> [1] MatrixGenerics_1.10.0 Biobase_2.58.0
#> [3] httr_1.4.4 sass_0.4.2
#> [5] foreach_1.5.2 bit64_4.0.5
#> [7] jsonlite_1.8.3 bslib_0.4.0
#> [9] assertthat_0.2.1 highr_0.9
#> [11] BiocManager_1.30.19 BiocFileCache_2.6.0
#> [13] blob_1.2.3 BSgenome_1.66.0
#> [15] GenomeInfoDbData_1.2.9 Rsamtools_2.14.0
#> [17] yaml_2.3.6 progress_1.2.2
#> [19] pillar_1.8.1 RSQLite_2.2.18
#> [21] lattice_0.20-45 glue_1.6.2
#> [23] digest_0.6.30 XVector_0.38.0
#> [25] htmltools_0.5.3 Matrix_1.5-1
#> [27] plyr_1.8.7 XML_3.99-0.12
#> [29] pkgconfig_2.0.3 biomaRt_2.54.0
#> [31] bookdown_0.29 zlibbioc_1.44.0
#> [33] BiocParallel_1.32.0 tibble_3.1.8
#> [35] KEGGREST_1.38.0 generics_0.1.3
#> [37] ellipsis_0.3.2 cachem_1.0.6
#> [39] SummarizedExperiment_1.28.0 GenomicFeatures_1.50.0
#> [41] cli_3.4.1 magrittr_2.0.3
#> [43] crayon_1.5.2 memoise_2.0.1
#> [45] evaluate_0.17 fansi_1.0.3
#> [47] doParallel_1.0.17 xml2_1.3.3
#> [49] tools_4.2.1 prettyunits_1.1.1
#> [51] hms_1.1.2 BiocIO_1.8.0
#> [53] lifecycle_1.0.3 matrixStats_0.62.0
#> [55] stringr_1.4.1 DelayedArray_0.24.0
#> [57] AnnotationDbi_1.60.0 Biostrings_2.66.0
#> [59] compiler_4.2.1 jquerylib_0.1.4
#> [61] rlang_1.0.6 grid_4.2.1
#> [63] RCurl_1.98-1.9 iterators_1.0.14
#> [65] rappdirs_0.3.3 VariantAnnotation_1.44.0
#> [67] rjson_0.2.21 bitops_1.0-7
#> [69] rmarkdown_2.17 restfulr_0.0.15
#> [71] codetools_0.2-18 curl_4.3.3
#> [73] DBI_1.1.3 R6_2.5.1
#> [75] GenomicAlignments_1.34.0 dplyr_1.0.10
#> [77] knitr_1.40 rtracklayer_1.58.0
#> [79] utf8_1.2.2 fastmap_1.1.0
#> [81] bit_4.0.4 filelock_1.0.2
#> [83] stringi_1.7.8 parallel_4.2.1
#> [85] Rcpp_1.0.9 vctrs_0.5.0
#> [87] png_0.1-7 tidyselect_1.2.0
#> [89] dbplyr_2.2.1 xfun_0.34