hpages
commented
2 years ago Hi,
A genomic range can be out-of-bound at the beginning or at the end of the chromosome, or at both ends. If the GRanges object "knows" the length of the chromosome, and if the chromosome is not circular, a warning gets issued when the object is transformed into an object that contains out-of-bound ranges:
library(GenomicRanges)
gr <- GRanges("chr1:55-80", seqlengths=c(chr1=100))
seqinfo(gr)
# Seqinfo object with 1 sequence from an unspecified genome:
# seqnames seqlengths isCircular genome
# chr1 100 NA <NA>With shift():
gr2 <- shift(gr, 30)
# Warning messages:
# 1: In valid.GenomicRanges.seqinfo(x, suggest.trim = TRUE) :
# GRanges object contains 1 out-of-bound range located on sequence chr1.
# Note that ranges located on a sequence whose length is unknown (NA) or
# on a circular sequence are not considered out-of-bound (use
# seqlengths() and isCircular() to get the lengths and circularity flags
# of the underlying sequences). You can use trim() to trim these ranges.
# See ?`trim,GenomicRanges-method` for more information.
# 2: In valid.GenomicRanges.seqinfo(x, suggest.trim = TRUE) :
# GRanges object contains 1 out-of-bound range located on sequence chr1.
# Note that ranges located on a sequence whose length is unknown (NA) or
# on a circular sequence are not considered out-of-bound (use
# seqlengths() and isCircular() to get the lengths and circularity flags
# of the underlying sequences). You can use trim() to trim these ranges.
# See ?`trim,GenomicRanges-method` for more information.
gr2
# GRanges object with 1 range and 0 metadata columns:
# seqnames ranges strand
# <Rle> <IRanges> <Rle>
# [1] chr1 85-110 *
# -------
# seqinfo: 1 sequence from an unspecified genomeWith promoters():
prom1 <- promoters(gr, upstream=60, downstream=0)
# Warning message:
# In valid.GenomicRanges.seqinfo(x, suggest.trim = TRUE) :
# GRanges object contains 1 out-of-bound range located on sequence chr1.
# Note that ranges located on a sequence whose length is unknown (NA) or
# on a circular sequence are not considered out-of-bound (use
# seqlengths() and isCircular() to get the lengths and circularity flags
# of the underlying sequences). You can use trim() to trim these ranges.
# See ?`trim,GenomicRanges-method` for more information.
prom2 <- promoters(gr, upstream=0, downstream=50)
# Warning message:
# In valid.GenomicRanges.seqinfo(x, suggest.trim = TRUE) :
# GRanges object contains 1 out-of-bound range located on sequence chr1.
# Note that ranges located on a sequence whose length is unknown (NA) or
# on a circular sequence are not considered out-of-bound (use
# seqlengths() and isCircular() to get the lengths and circularity flags
# of the underlying sequences). You can use trim() to trim these ranges.
# See ?`trim,GenomicRanges-method` for more information.Note that these warning messages actually come from the validity method for GRanges objects:
validObject(gr2)
# [1] TRUE
# Warning message:
# In valid.GenomicRanges.seqinfo(x, suggest.trim = TRUE) :
# GRanges object contains 1 out-of-bound range located on sequence chr1.
# Note that ranges located on a sequence whose length is unknown (NA) or
# on a circular sequence are not considered out-of-bound (use
# seqlengths() and isCircular() to get the lengths and circularity flags
# of the underlying sequences). You can use trim() to trim these ranges.
# See ?`trim,GenomicRanges-method` for more information.
validObject(prom1)
# [1] TRUE
# Warning message:
# In valid.GenomicRanges.seqinfo(x, suggest.trim = TRUE) :
# GRanges object contains 1 out-of-bound range located on sequence chr1.
# Note that ranges located on a sequence whose length is unknown (NA) or
# on a circular sequence are not considered out-of-bound (use
# seqlengths() and isCircular() to get the lengths and circularity flags
# of the underlying sequences). You can use trim() to trim these ranges.
# See ?`trim,GenomicRanges-method` for more information.
validObject(prom2)
# [1] TRUE
# Warning message:
# In valid.GenomicRanges.seqinfo(x, suggest.trim = TRUE) :
# GRanges object contains 1 out-of-bound range located on sequence chr1.
# Note that ranges located on a sequence whose length is unknown (NA) or
# on a circular sequence are not considered out-of-bound (use
# seqlengths() and isCircular() to get the lengths and circularity flags
# of the underlying sequences). You can use trim() to trim these ranges.
# See ?`trim,GenomicRanges-method` for more information.Also note that these objects can easily be trimmed with trim() (as explained in the above warning messages):
trim(gr2)
# GRanges object with 1 range and 0 metadata columns:
# seqnames ranges strand
# <Rle> <IRanges> <Rle>
# [1] chr1 85-100 *
# -------
# seqinfo: 1 sequence from an unspecified genome
trim(prom1)
GRanges object with 1 range and 0 metadata columns:
# seqnames ranges strand
# <Rle> <IRanges> <Rle>
# [1] chr1 1-54 *
# -------
# seqinfo: 1 sequence from an unspecified genome
trim(prom2)
# GRanges object with 1 range and 0 metadata columns:
# seqnames ranges strand
# <Rle> <IRanges> <Rle>
# [1] chr1 55-100 *
# -------
# seqinfo: 1 sequence from an unspecified genomeThis is quickly mentioned in ?GenomicRanges::trim.
The reason for the current design is that negative positions are not necessarily a problem. They are actually meaningful on a circular chromosome like chrM:
library(BSgenome.Hsapiens.UCSC.hg38)
getSeq(BSgenome.Hsapiens.UCSC.hg38, GRanges("chr1", IRanges(-4, 10)))
# Error in loadFUN(x, seqname, ranges) :
# trying to load regions beyond the boundaries of non-circular sequence
# "chr1"
chrM <- BSgenome.Hsapiens.UCSC.hg38$chrM
chrM
# 16569-letter DNAString object
# seq: GATCACAGGTCTATCACCCTATTAACCACTCACGGG...AGCCCACACGTTCCCCTTAAATAAGACATCACGATG
getSeq(BSgenome.Hsapiens.UCSC.hg38, GRanges("chrM", IRanges(1, 10)))
# DNAStringSet object of length 1:
# width seq
# [1] 10 GATCACAGGT
getSeq(BSgenome.Hsapiens.UCSC.hg38, GRanges("chrM", IRanges(-4, 10)))
# DNAStringSet object of length 1:
# width seq
# [1] 15 CGATGGATCACAGGTSo it's not clear that they should be automatically and systematically "fixed". Furthermore, doing so would make some operations hard to revert. For example shift(gr, -30) would not longer be guaranteed to revert shift(gr, 30).
This is why handling out-of-bound ranges is the responsibility of the user.
I understand that this is maybe harder to do in the context of getPromoterSeq but you're not providing a reproducible example so it's hard for me to further comment on that.
H.
cmatKhan
Should
promoters()be allowed to return negative indicies?When trying to use
getPromoterSeq(), I received the following error:Upon inspecting the promoter regions with the same setting ,
upstream = 400, as I had used ingetPromoterSeq, I discovered that there are two genes which, with this promoter definition, return negative values due to the fact that they are annotated at the very beginning of the chromosome such that 400bp upstream causes the promoter range to be returned as a negative value:I confirmed that this was, in fact, the cause of the error in
getPromoterSeq()by settingupstreamto 20 or less, so that no negative values are returned, andgetPromoterSeq()returned without error.Possibly in the
promotersfunction, a warning should be issued? Alternatively, if the definition of the promoter places the upstream/downstream indicies out of range of the chromosome, maybe setting to the chromosome limit (0 or length(chromosome)) would be appropriate? Maybe this should be an argument, something like,handle_out_of_bound_indicies = TRUEby default?