GRange: Average Metadata Value from 3 Replicates

Question

[vanbelj]

I have CPM-normalized genome coverage data from three replicate experiments in the form of bigwig files. I read these into R using rtracklayer and now have 3 GRange objects.

For each basepair, I want to calculate the mean score of the three replicates and save this as a new GRange object. I'll the GRange with averaged values as input for the EnrichedHeatmap package.

Is there a way to do this using GenomicRanges? If not, how would you go about getting average coverage data from 3 biological replicates?

> Rap1_Rep1
GRanges object with 1781467 ranges and 1 metadata column:
            seqnames      ranges strand |     score
               <Rle>   <IRanges>  <Rle> | <numeric>
        [1]     chrI       1-803      * |   0.00000
        [2]     chrI     804-805      * |   2.01756
        ...      ...         ...    ... .       ...

> Rep2
GRanges object with 1627264 ranges and 1 metadata column:
            seqnames      ranges strand |     score
               <Rle>   <IRanges>  <Rle> | <numeric>
        [1]     chrI       1-803      * |   0.00000
        [2]     chrI     804-809      * |   2.02448
        ...      ...         ...    ... .       ...

> Rep3
GRanges object with 1537131 ranges and 1 metadata column:
            seqnames      ranges strand |     score
               <Rle>   <IRanges>  <Rle> | <numeric>
        [1]     chrI       1-803      * |   0.00000
        [2]     chrI     804-813      * |   3.61258
        ...      ...         ...    ... .       ...

Answer 1

[vanbelj]

The process I developed to do this in R was extremely slow. I ended up calculating the mean of the bigwig files before importing into R by following ATpoint's response here.

Answer 2

[Jeremy Leipzig]

some of your ranges aren't comparable but you can disjoin them and hopefully the metadata will survive

grl <- GRangesList(rep1,rep2,rep3)
disjoin(grl)

[vanbelj]

Hi Jeremy,

This is a result of the compression done by bigwig files. I can extract the coverage data for each basepair to a list with the code below. I can then performe the averages inside a dataframe but I have no clue how to get this information back into a GRange.

Rep1_cov <- rep(Rep1$score, width(Rep1))
Rep2_cov <- rep(Rep2$score, width(Rep2))
Rep3_cov <- rep(Rep3$score, width(Rep3))

> length(Rep1_cov)
[1] 12157105
> length(Rep2_cov)
[1] 12157105
> length(Rep3_cov)
[1] 12157105