I have microarray analysis problem that I need to solve in R. Mostly it is data and genomic ranges manipulation in R, but my R skills are not so good, hence I am seeking help from the community.
What I want:
To merge (bin) microarray probes into windows and calculate statistics for every window.
Why I want this:
To do cluster analysis for bined probes.
Example of data:
> probeCoordintes
Chr probeStart probeEnd
P1 chr7 4655 4680
P2 chr7 4691 4716
P3 chr7 4724 4749
P3 chr7 4757 4782
P4 chr7 4787 4812
P5 chr7 4824 4849
> Norm.Probes
Case1 Case2 Case3
P1 12.57384 12.17727 12.40320
P2 13.42015 12.73005 12.71374
P3 11.77102 11.41472 11.46454
P4 11.74266 11.39161 11.52613
P5 11.62942 11.29788 11.59295
Workflow of what I need:
- Merge probes into windows of x bp size (e.g, 1kb) according coordinates
- For every Case, For every Window calcute statistics (e.g, mean)
- Perfect result would look like this:
windows
W1 chr1:0-1000 W2 chr1:1000-2000 ...
values_means
Case1 Case2 Case3
W1 11 12 13
W2 10 10 11
...
I know I am asking much without any code of my own, but hope that someone will help me with this.
1 answer
First, you should get your data in an appropriate data structure, such as a GRanges:
library(GenomicRanges)
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
probes <- with(probeCoordinates,
GRanges(Chr, IRanges(probeStart, probeEnd),
seqlengths=seqlengths(TxDb.Hsapiens.UCSC.hg19.knownGene)))
Next, you should form the GRanges representing the windows:
tiles <- unlist(tileGenome(seqinfo(probes), tilewidth=1000L))
Then, find the overlaps between the probes and windows and take the mean according to the overlap grouping. This is a little bit hairy but it should be efficient.
hits <- findOverlaps(probes, tiles)
s <- rowsum(as.matrix(Norm.Probes[queryHits(hits),]), subjectHits(hits))
counts <- countSubjectHits(hits)
presentWindows <- counts > 0L
tiles$means <- matrix(0L, ncol=ncol(Norm.Probes), nrow=length(tiles))
tiles$means[presentWindows,] <- s / counts[presentWindows]
Does this look like what you want?
subset(tiles, rowSums(means) > 0)
GRanges with 2 ranges and 1 metadata column:
seqnames ranges strand | means
<Rle> <IRanges> <Rle> | <matrix>
[1] chr7 [3685, 4684] * | 12.57384 12.17727 12.4032
[2] chr7 [4685, 5684] * | 12.1408125 11.708565 11.82434
---
seqlengths:
chr1 chr2 ... chrUn_gl000249
249250621 243199373 ... 38502
•
written
by
Michael Lawrence •
Log in to answer this question.
Have you looked into the summarizeOverlaps function of Bioconductors GenomicRanges library?
<deleted>