Hello everyone,
I am using the GenomicRanges module in R to output the maximum and minimum of my genomic range object.
First, the GRanges object has been designed with a file columns. Note that metadata have been added to the object using "gc=".
> gr <- GRanges(seqnames=file$V1, > ranges=IRanges(start=as.numeric(file$V7), end=as.numeric(file$V8)), > gc=file$V2)
Unfortunately, using the range() function, the metadata column are lost.
gr_ran <- range(gr, ignore.strand = TRUE)
#> gr_ran
#GRanges object with 1024 ranges and 0 metadata columns:
# seqnames ranges strand
# <Rle> <IRanges> <Rle>
# [1] Scaffolds_0_pilon 327519-354933 *
# [2] Scaffolds_7_pilon 228693-230062 *
Would anyone know how to keep these metadata column?
You cannot keep the metadata in the way you running your command.
‘range’ returns an object of the same type as ‘x’ containing range bounds for each distinct (seqname, strand) pairing. The names (‘names(x)’) and the metadata columns in ‘x’ are dropped.
The range call will collapse your ranges by strand+chromosome and thus there is no obvious or straightforward way to transfer the meta.data.
For example, if you have 2 regions to collapse each with their own metadata which value do you take? And then what if there are 3 regions, etc...
If you are looking for the start and end of each region you can use start() and end(). [keep in mind strandedness].
If you want to do exactly what you are seem to be explaining I would suggest splitting the whole GRanges object by chromosome+strand followed by taking the mean of the gc column using lapply or dplyr or something of that nature.
Thank you very much for your answer benformatics. You have well understood my problem, and in my case the gc column is not numbers, so there is indeed no way to make a mean of it.
Log in to answer this question.
